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technique  for  determining  the  aerodynamic  parameters 
of  a  maneuvering  reentry  vehicle  (MaRV)  from  post-flight 
analys.'S  of  radar  tracking  data  was  developed.  The  techni¬ 
que  uses  a  sequential  weighted  least  squares  estimator  which 
steps  through  the  data  processing  batches  of  data  sequen¬ 
tially.  The  observations  consisted  of  angle  and  range  measure 
ments  from  a  precision  tracking  radar  located  in  the  vicinity 
of  the  impact  area.  A  six  dimension  estimator  consisting  of 
three  position  and  three  velocity  components  was  designed 
for  estimating  the  state  of  the  RV  during  free  flight.  Two 
reentry  estimators  were  developed.  One,  a  seven  dimension 
estimator  consisting  of  the  six  elements  of  the  free-f light 
estimator  plus  the  ballistic  coefficient,  was  designed  to 
estimate  the  state  during  ballistic  reentry.  The  second  is 
the  nine  element  MaRV  estimator  used  once  a  maneuver  is  be¬ 
gun.  An  algorithm  based  on  measurement  residual  monitoring 
was  dev'feloped  to  switch  adaptively  from  the  six-state  esti¬ 
mator  to  the  seven-state  estimator  and  then  to  the  nine-state 
estimator.  A  series  of  numerical  simulations  was  performed 
to  validate  the  technique  and  its  programming.  Monte  Carlo 
simulations  were  used  to  verify  the  accuracy  of  the  estima¬ 
tor  covariance  matrix. 


POST -FLIGHT  TRAJECTORY  RECONSTRUCTION  OF  A 
MANEUVERING  REENTRY  VEHICLE  FROM 
RADAR  MEASUREMENTS 

I  Introduction 

Objective  and  Scope 

The  objective  of  this  thesis  was  to  develop  a  technique 
for  the  efficient  and  accurate  post-mission  analysis  of  radar 
tracking  data  collected  on  the  reentry  vehicle  (RV)  of  a 
ballistic  missile.  The  RV  was  assumed  to  have  the  capability 
of  maneuvering  in  the  atmosphere  by  controlling  the  direction 
and  magnitude  of  its  lift  vector.  This  ignores  the  case  of  a 
maneuverable  reentry  vehicle  (MaRV)  which  might  be  designed 
to  maneuver  by  thrusting}  and  therefore,  could  conceivably 
maneuver  even  above  the  atmosphere.  In  general,  the  technique 
which  was  developed  and  which  is  described  in  this  report 
would  not  be  applicable  to  such  an  RV.  It  was  further  assumed 
that  a  single  radar  was  ideally  located  in  the  vicinity  of  the 
impact  area  and  that  the  radar  tracked  the  RV  from  near  hori¬ 
zon  break  (a/O®  elevation)  to  very  near  impact.  The  accura¬ 
cies  of  the  radar  tracking  data,  which  consists  of  time-tagged 
range,  azimuth,  and  elevation  measurements,  were  assumed  to 
be  consistent  with  those  of  a  state-of-the-art,  precision, 
phased-array  tracking  radar.  The  data  rates  were  also  con¬ 
sidered  to  be  within  the  capability  of  such  a  radar.  ATI  of 
the  data  which  were  used  for  this  investigation  were 
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simulated.  Table  I  lists  the  pertinent  information  relative 
to  the  simulated  data  where  the  noise  on  the  data  was  simu¬ 
lated  as  zero-mean  white  Gaussian  noise. 

Table  I-l 

Description  of  Simulated  Data  Used  in  Numerical  Examples 


Data 

Rates  t 

Range 

1  point  per  second  above 

Azimuth 

120  km  and  10  points  per 

Elevation 

second  below  120  km 

Data  1 

-  a  Errors 

Range 

5  m 

Azimuth 

0.016® 

Elevation 

0.020® 

A  further  assumption  is  that  all  parameters  other  than 
those  being  estimated  are  perfectly  Known.  This  is  equivalent 
to  saying  that  there  are  no  "Q-parameters",  where  Q-parameters 
are  defined  by  Day  (Ref  11)  and  others  as  those  parameters 
affecting  the  observations  but  which,  for  some  reason  or 
another,  are  not  or  cannot  be  estimated.  These  include, 
among  others,  station  location,  biases,  and  a  most  important 
one  for  the  technique  developed,  atmospheric  density.  This 


is  not  a  limitation  on  the  technique  itself  but  on  the  inter 
pretation  of  the  covariance  matrices  which  result. 


The  Research  Approach 


The  problem  of  real-time  tracKing  of  a  reentry  vehicle 
based  on  its  radar  measurements  has  received  considerable  at¬ 
tention  in  the  literature  over  the  past  2-3  decades.  It  has 
been  described  by  various  authors  (25)  and  (10)  as  a  highly 
complex  problem  in  nonlinear  filtering.  Since  1960,  following 
the  important  work  of  Kalman  and  Bucy  (Ref  20),  the  extended 
Kalman  filter  (EKF)  has  found  wide  application  for  both  the 
ballistic  reentry  vehicle  (BRV)  tracking  problem  (22),  (14) 
and  later  the  maneuvering  reentry  vehicle  tracking  problem 
(10).  This  is  owing  to  the  real  time  nature  of  the  problem 
where  measurements  need  to  be  processed  as  they  become  avail¬ 
able  so  as  to  "point**  the  radar  for  the  next  measurement. 

Here  the  primary  emphasis  is  on  maintaining  track  on  the  target. 

The  post-flight  analysis  problem  is  very  similar  to  the 

tracking  problem,  but  it  differs  in  two  important  aspects. 

The  similarity  is  that  the  true  equations  of  motion  governing 

the  flight  of  the  KV  being  tracked  are  the  same  in  both  cases. 

The  differences  involve  1),  the  time  constraint  and  2),  the 

emphasis.  In  the  post-flight  analysis,  there  is  virtually  no 

time  constraint  and  the  emphasis  is  on  accuracy.  That  is  to 

say  that  we  are  interested  in  determining  the  RV  state  varia- 

* 

bles  at  some  epoch  time  (or  possibly  a  series  of  epoch  times) 
with  as  much  accuracy  as  the  data  will  permit.  It  should, 
after  all,  be  the  data  which  limit  the  accuracy  of  our  esti¬ 
mates  and  not  our  model.  In  the  tracking  problem,  out  of 
necessity  imposed  by  the  time  constraint,  many  simplifying 

Defined  as  the  time  at  which  state  variables  are  estimated. 
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assumptions  about  the  dynamics  model  are  made,  o«g.,  a  flat, 
nonrotating  earth  may  be  assumed.  In  the  post-mission 
emalysis  problem,  even  a  rotating,  spherical  earth  is  not  of 
sufficient  accuracy.  Because  of  this  emphasis  on  accuracy, 
much  attention  to  detail  is  given  in  the  later  sections  of 

this  report  which  deal  with  modeling  the  flight  of  the  MaRV. 

■ »  »  . 
It  should  be  pointed  out  that  the  term  state  of  the  RV  as  used 

above  implies  the  RV  state  vector  which  consists  of  position, 
velocity,  and  zero,  one,  or  three  aerodynamic  parameters  de¬ 
pending  respectively  on  whether  the  RV  is  in  free  flight, 
ballistic  flight,  or  maneuvering  flight.  Quite  often,  in 
missile  analysis,  the  aerodynamic  parameters  are  the  ones  ve 
are  most  interested  in  estimcu.'ng  in  post-mission  trajectory 
reconstruction . 

Another  advantage  we  have  in  the  post-flight  analysis 
problem  over  the  real-time  tracking  problem  is  that  we  have 
all  the  data  at  hand  to  treat  simultaneously,  if  ve  desirei 
or  one  point  at  a  time,  if  we  desire.  This  gives  us  great 
flexibility  in  designing  the  estimator  for  the  post-flight 
analysis  problem.  It  also  forces  us  to  make  a  choice.  Chang 
et  al.,  (Ref  8)  have  applied  a  fixed  interval  smoother,  vhich 
consists  of  combining  the  outputs  of  a  forward  pass  of  an  ex¬ 
tended  Kalman  filter  with  a  backward  pass  of  the  extended 
Kalman  filter,  to  solve  the  problem.  Gauss  invented  the 
method  of  weighted  least  squares  (WiS)  in  1795  to  solve  an 
astrodynamical  problem  by  treating  all  of  the  available  obser¬ 
vations  simultaneously.  The  technique,  which  was  developed 
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for  thif  roioarch  i«  a  vary  flaxibla  algorithm  vhich  ia  basad 
on  tha  mathod  of  WLS.  Tha  uaar«  by  apecifying  tha  batch  aiza 
(if a. I  tha  numbar  of  obaarvationa  to  procaaa  in  a  givan  batch )« 
can  oparata  tha  program  in  aithar  a  batch  VfLS  moda  vhara 
all  data  ara  traatad  aimultanaouaiyi  or  aa  an  EKP  traating 
aach  data  point  (ranga,  azimuth,  and  alavation  at  a  aingla 
tima)  aaparatalyj  or  any  placa  in  batvaon  aa  a  aaquantial  WL8 
aatimator. 

Tha  algorithm  vhich  vaa  davalopad  and  programmad  for 
thia  thaaia  raaaarch  ia  praaantad  in  datail  in  tha  ramaining 
chaptara  of  thia  raport«  it  haa  two  noval  faaturaa.  Ona  ia 
tha  way  in  vhich  tha  trajactory  partiala,  vhid)  ara  raquirad 
in  tha  vaightad-iaaat-aquaraa-diffarantial-oorraction  procaaa, 
ara  calculated #  Tha  other  ia  a  taohniqua  for  adaptively 
avitching  from  tha  aix-alamant  atata  vector  uaad  for  fraa- 
f light  aatimation  to  tha  aavan-alamant  atata  vector  uaad  for 
balliatic-raantry-flight  and  than  from  tha  aavan-alamant  atata 
vector  to  tha  nina-alamant  MaAV  atata  vector # 

Tha  atandard  technique  for  calculating  trajactory  par- 
tiala  in  both  aatallita  and  miaaila  trajectory  raconatruction 
ia  to  integrate  numerically  a  eat  of  aacond  order,  ordinary 
differential  aquationa  called  variational  aquationa.  8aa  for 
example  Rafa  (1),  (27)  and  (34),  Tha  varlationaJ  aquationa 
have  aa  their  dependant  variablaa  tha  partial  darlvativaa  of 
tha  total  acceleration  with  raapact  to  tha  paramatara  to  be 
aotimatadf  Thera  ara  at  laaat  (fraa-f light)  10  of  thaaa  aacond 


order  (corresponding  to  36  first  order)  differential  equations 
to  be  integrated.  They  are  typically  passed  to  the  same  in¬ 
tegration  routine  which  Integrates  the  equations  of  motion. 

This  represents  a  significant  computational  burden  on  a  tech- 
niquSf  whether  batch  or  sequential,  which  is  at  best  slow) 
however,  because  the  WLS  differential  correction  technique  is 
inherently  iterative,  the  accuracies  of  the  partials  need  not 
be  as  great  as  those  of  the  state  variable  themselves .  The 
validity  of  this  claim  is  demonstrated  in  Chapter  IV  of  this 
report  where  it  is  shown  that  the  batch  least  squares  estima¬ 
tor  converges  in  two  or  three  iterations  with  good  initial, 
estimates  of  the  state  variables. 

As  an  example,  a  spherical  harmonic  series  of  zonal 
harmonics,  representing  the  potential  of  an  oblote  spheroid, 
is  used  for  accurately  modeling  the  missile’s  acceleration  due 
to  gravity)  whereas,  the  predominate  term  of  this  series,  that 
due  to  a  spherical  earth,  is  quite  adequate  for  partials  cal¬ 
culations.  As  s  result  of  this  less  demanding  requirement  on 
the  partial^  accuracies,  a  very  simple  Taylor’s  series  expan¬ 
sion  is  used  to  equate  the  state,  at  some  time  with 

the  state,  2^^  at  some  earlier  time  (i.e.,  numerically  in- 
tegrate  This  is  then  used  in  calculating  the  state  tran¬ 

sition  matrix,  j5(ti^j,  t^)  which  is  central  to  the  trajectory 
partials*  calculations  as  well  as  to  propagating  the  state 
covariance  matrix  in  time.  The  approximations  mentioned  above, 
as  well  as  others  used  in  the  partials  calculations,  are  given 
in  complete  detail  in  Appendix  fi  along  with  all  the  partials 


used  in  the  development  of  the  state  transition  matrix.  It 
should  be  emphasized  that  these  approximations  are  only  used 
in  partials  calculations. 

The  other  unique  feature  of  the  technique  which  was 
developed  is  a  method  for  adaptively  sensing  when  a  higher 
order  estimator  is  required  and  then  switching  to  the  higher 
order  estimator.  As  was  mentioned  earlier »  it  was  assumed 
that  the  data  collection  began  as  the  RV  broke  the  radar's 
horizon  mecUiing  that  several  mim  i  <• free-flight  data  were 
collected.  Free-flight  is  that  phase:  of  the  trajectory  after 
booster  burnout  and  prior  to  reentry  when  the  RV  is  acceler¬ 
ated  only  by  the  force  of  the  earth's  gravity.  During  this 
phase  of  flight,  the  RV  state  vector  is  given  as 

X  »  (x,y,z4,y,z)'^ 

Where  X,  Y,  and  Z  are  coordinates  of  the  RV's  position  vector 
and  X,  Y,  and  Z  are  coordinates  of  the  RV’s  velocity  in  some 
arbitrarily  chosen  cartesian  coordinate  system.  In  this  con¬ 
text,  the  state  vector  is  a  vector  of  those  parameters  we 
wish  to  estimate  at  some  given  time  called  the  epoch  time 
such  that  given  the  state  vector  at  epoch  and  the  equation  of 
motion  for  each  state  variable,  we  can  predict  the  state  at 
any  other  time  within  the  interval  during  which  the  equations 
are  valid.  Alter  reentering  the  earth's  atmosphere,  the  RV 
begins  to  experience  an  additional  force,  that  due  to  atmos¬ 
pheric  drag.  During  this  ballistic-reentry-phase*  the  RV's 

♦Even  a  MaRV  initially  reenter  along  a  ballistic  tra¬ 
jectory.  The  MaRV  state  vector  is  given  in  later  sections 
of  this  report. 
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state  vector  is  given  as 


X  «  (X.Y,Z,X.Y,Z,3)'^ 

where  3  is  a  parameter  associated  with  the  drag  called  the 
ballistic  coefficient.  During  free-flight,  the  ballistic 
coefficient  is  completely  unobservable  (meaning  it  can  not  be 
estimated  from  the  data) .  The  transition  from  free-flight  to 
reentry,  where  P  is  observable,  is  a  gradual  one  with  no 
clear  cut  demarcation}  therefore,  if  one  is  using  a  sequential 
estimator  to  step  through  the  data,  there  must  be  a  way  to 
determine  when  this  transition  from  free-flight  to  ballistic- 
reentry-flight  has  occurred.  In  chapter  III  of  this  report, 
a  method  based  on  monitoring  the  root -mean-square  (RMS)  error 
of  the  data  residuals  is  presented  for  this  purpose. 

Report  Organization 

The  remaining  sections  of  this  report  describe  in  detail 
the  necessary  calculations  for  the  post-flight  analysis  of 
radar  tracking  data  collected  on  a  maneuvering  reentry  vehi¬ 
cle.  Chapter  II  presents  the  equations  necessary  for  accu¬ 
rately  modeling  the  flight  of  a  MaRV.  These  include  the  equa¬ 
tions  of  motion,  various  transformations  need,  and  models  of 
the  earth's  atmosphere  and  the  RV's  drag.  In  Chapter  III, 
the  algorithm  for  estimating  the  trajectory  and  aerodynamic 
parameters  from  the  noise  corrupted  radar  measurements  is 
developed.  This  starts  with  the  basic  weighted  least  squares 
(WLS)  estimator  for  nonlinear  equations  and  shows  how  this  is 
related  to  the  minimum  variance  (MV)  and  maximum  likelihood 
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(ML)  estimators .  A  sequential  estimator  with  and  without 
a  priori  covariance  on  the  initial  conditions  is  presented. 

A  technique  for  adaptively  switching  from  the  free-f light 
estimator  to  the  BRV  estimator,  to  the  MaRV  estimator  is 
given,  and  a  method  for  determining  close  initial  estimates 
is  suggested.  Chapter  IV  discusses  the  numerical  simulations 
that  went  into  developing  and  testing  the  computer  program 
which  was  developed.  Two  numerical  examples  are  presented 
in  detail  which  demonstrate  the  operation  of  the  program 
using  simulated  data.  Chapter  V  presents  a  brief  summary 
and  suggests  some  areas  where  further  study  is  required.  Two 
appendices  are  also  included.  Appendix  A  describes  the 
oblate  spheroid  which  is  used  as  the  model  of  the  geometric 
shape  of  the  earth.  Appendix  B  defines  the  state  transition 
matrix  and  presents  the  equations  used  in  deriving  the  par¬ 
tial  derivatives  of  the  state  transition  matrix,  and  then 
lists  all  the  partials  thus  derived. 
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II  Trajectory  Generation  Equations  and  Methods 


This  chapter  presents  the  equations  necessary  for  ac¬ 
curately  modeling  the  flight  of  a  ballistic  missile's  MaPV, 

These  include  the  differential  equations  of  motion  and  a 
numerical  integration  scheme  for  propagating  the  MaRV's  state 
vector  in  time.  Also  included  are  the  equations  for  trans¬ 
forming  the  position  of  the  MaRV  in  the  computational  co¬ 
ordinate  system,  the  one  in  which  the  equations  of  motion  are 
expressed  and  integrated,  to  radar  range,  azimuth,  and  eleva¬ 
tion  for  a  given  radar  station.  This  transformation  will  be 
necessary  for  the  least  squares-dif ferential  correction 
algorithm  presented  in  chapter  III  which  compares  the  observed 
radar  measurements  to  the  calculated  radar  measurements  based 
on  an  assumed  or  estimated  state  vector  for  the  MaRV.  Also 
included  in  this  chapter  are  models  of  the  eartHs  atmosphere 
and  the  drag  model  for  a  class  of  sphere-cone  shaped  RV's. 

Coordinate  Systems  and  Transformations 

Several  different  coordinate  systems  are  useful  for 
expressing  the  trajectory  and  radar  related  equations.  The 
most  important  ones  are  described  in  this  section.  The  com¬ 
putational  coordinate  system,  the  one  in  which  the  equations 
of  motion  are  expressed  and  integrated,  is  the  earth-centered- 
inertial  (ECI)  cartesian  coordinate  system  depicted  in  Figure  II- 
1.  The  origin  of  this  coordinate  system  is  at  the  center  of 
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the  ellipsoidal  earth  model *j  the  X-Y  plane  is  the  equatorial 
plane  of  the  earthi  and  the  Z-axis  is  coincident  with  the 
earth's  axis  of  rotation.  This  coordinate  system  is  fixed  in 
inertial  space  with  the  X-axis  passing  through  the  Greenwich 
meridian  at  some  arbitrarily  chosen  initialization  time,  tj. 
For  this  analysis,  tj  is  taken  as  the  time  of  the  first  radar 
observation. 

Figure  II -1  also  shows  the  radar-RV  geometry.  The 
radar  location  is  typically  expressed  in  geographic  coordi- 

iff 

nates -geodetic  latitude  (^  ),  longitude  (^^),  and  height  (h) 
above  the  earth's  surface.  These  coordinates  are  also  con¬ 
venient  for  displaying  and  plotting  the  RV's  position  as  a 
function  of  time.  Figure  II-2  identifies  the  geographic  co¬ 
ordinates  of  the  radar  station. 

Two  other  coordinate  systems  which  are  useful  in  de¬ 
scribing  the  trajectory  estimation  problem  are  two  radar  co¬ 
ordinate  systems  which  are  shown  in  Figure  II-3.  One  is  a 
left-handed  cartesian  coordinate  system  (x,y,z)  called  the 
north-east-up  (NEU)  system  and  the  other  is  a  spherical  co¬ 
ordinate  system  (R,A;5,E).  The  origin  of  these  two  coordinate 
systems  is  the  radar.  The  z-axis  is  along  the  local  geodetic 
vertical  and  azimuth  an  elevation  are  measured  as  shown  in 
the  figure. 

The  final  coordinate  system  to  be  described  in  this 
section  is  the  earth-centered-rotating  (ECR)  system.  It 


.V 


*See  Appendix  A  for  a  description  of  the  reference  el¬ 
lipsoid  used  to  model  the  earth. 
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FIG  II-3  GEOCENTRIC  COORDINATES  OF  THE  RADAR  STATION 


13 


differs  from  the  ECI  coordinate  system  only  in  that  the  X-¬ 
axis  alvays  points  out  the  Greenwich  meridian.  The  ECI  and 
ECR  coordinate  are  related  by  the  following  transformation  at 
any  given  time  t^ . 


'x,' 

y. 

ccs(^ifC)  o 

y 

z,. 

o  o  / 

2 

— 

(II-l) 


where 


At  =  t.-4 


n  is  the  earth's  rotation  rate  and 

Xp,  Yp,  Zp  are  coordinates  in  the  ECR  reference  frame 

The  geographic  coordinates  h)  are  related  to 

the  ECR  coordinates  by  the  following  equations. 


cas  ^  A 


(II-2) 


ihi)  (II-3) 

h  -  (II-5) 

=  swXz/r)  (II-6) 

A’  OiA'XYr  / Xk) 
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and 


s  (II-8) 

€  '  (II-9) 

where 

^  is  the  geocentric  latitude 

is  the  local  radius  of  the  earth 
a  is  the  semi -major  axis  of  the  ellipsoidal  earth  model 
b  is  the  semi-minor  axis  of  the  ellipsoidal  earth  model 
and 

e  is  the  eccentricity  of  the  ellipsoidal  earth  model. 
Equations  (II-8)  and  (II-9)  are  derived  in  Appendix  A,  Appen 
dix  A,  which  describes  in  detail  the  reference  ellipsoid  used 
to  model  the  shape  of  the  earth,  also  gives  the  numerical 
values  of  the  constants  a,  b,  and  e. 


The  final  transformation  which  will  be  necessary  to 
transform  the  ECI  coordinates  of  the  RV  to  radar  observation 
coordinates  (R,  Az,  and  E)  is  the  following 


f 

- 

-5/A' A, 

005 

o 

%-x 

cos 

ccs^i  s/^'Xi 

z-z. 

where 

^  , 
s' 

and  h„  are  the 
s 

geographic  coordinates 

of  the 

radar  station  (Fig  I I -2  )  and  X  ,  y  and  Z  are  the 

w  s  s 
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coordinates  of  the  radar  station  in  t>ie  ECR  reference  frame 
given  by 


^5  ^ 

(II-ll) 

CA-^^i)caS^s 

(11-12) 

^9^  C^s 

(11-13) 

where  r  is  the  radius  of  the  earth  at  the  radar  station  and 
es 

is  the  geocentric  latitude  of  the  radar  station.  and 


* 

r^g  are  determined  from^^  using 


Then 

Ai  *  //x) 

/=  wi:  >■//?) 


(II-8)  and  (II-9). 

(11-14) 

(11-15) 

(11-16) 


Free-Fliqht  Equations  of  Motion 

For  a  reentry  vehicle,  unpowered  and  still  above  the 
earth's  atmosphere,  the  only  significant  force  acting  on  it 
is  that  due  to  the  earth's  gravity.  The  acceleration  due  to 
gravity  is  typically  derived  from  the  gravity  potential. 
Wolaver  (Ref  38)  shows  that  the  earth's  gravity  potential, 

U  is  a  solution  to  Laplaces  equation 


(11-17) 


tu  .  tu 

ay*- 
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The  roost  general  solution  to  Eq  (11-17)  for  a  unit  mass  ex¬ 
terior  to  the  earth  is  given  as  an  infinite  spherical  harmonic 
series  which  is  a  function  of  r,  the  distance  from  the  earth's 
center; the  geocentric  latitude;  and  X,  the  longitude 

V -It (K-18) 

where, 

is  the  earth's  mass  times  the  universal  gravity  con¬ 
stant 

a  is  the  mean  equitorial  radius 

M 

P^(sin<y),  Pj^  (sin^  are  the  Legendre  polynomials  and 

M  M 

^  are  constants  determined  from  satellite  data. 

Typically,  some  truncated  form  of  Eq  (11-18)  is  used  for  sat¬ 
ellite  and  missile  trajectory  simulation  and  analysis  depending 
on  the  accuracy  which  is  required,  A  particularly  convenient 
form  is  one  which  assumes  the  dependence  on  longitude  to  be 
negligible.  This  truncated  expression  (Ref  5) 

is  called  the  zonal  potential.  The  first  four  terms  of  the 
zonal  potential  were  used  in  this  investigation  to  model  the 


18 


earth's  gravity.  Values  of  the  constants  used  and  the  first 
three  Legendre  polynomials  are  given  belov 


(11-20) 

(11-21) 

(11-22) 

d  •  rm 

(11-23) 

/OSZ.6/6 

(11-24) 

55 

(XI-25) 

(11-26) 

^  *  596600.  B  kn^* /it;.*' 

(11-27) 

The  zonal  potential  is  considered  sufficiently  accurate  for 
the  application  being  considered  hers  and  the  advantage  is 
that  by  setting 

5/A^f  *  ^/r  (11-28) 

it  can  be  vritten  completely  in  terms  of  the  computational 
coordinate  system.  This  eliminates  the  need  to  transform  from 
a  non-inertial  system  to  the  inertial  system  at  each  integra¬ 
tion  step  as  would  be  required  if  the  timo  dependent  terms 


v«r«  ratainttd*  Tha  diaadvantaga  of  using  this  potential  func- 

2 

tion  ia  that  tha  J2  tarm  is  naglactad.  The  magnitude  of  the 
2 

3 2  tarm>  vhich  represents  the  Earth'd  equatorial  ellipticity, 
has  tha  same  order  as  tha  Jj  zonal  harmonic  and  perturbs 
motion  as  much  as  any  zonal  harmonic  except  J2*  Thus  for 
other  applications,  tha  potential  raprasantad  by  Eq  (11-19) 
may  not  give  sufficient  accuracy. 

The  acceleration  due  to  gravity  is  given  as  tha  negative 
of  tha  gradient  of  tha  potential  function. 

(11-29) 

In  tha  ECl  rafaranca  frame,  tha  components  of  tha  accelera¬ 
tion  due  to  gravity  become 


(11-30) 


(11-31) 


V Z(0f(:  15 * 70^  -nf)J  (11-32) 

The  RV  stats  vector  during  free  flight  than  contains  tha  six 
elements  of  position  and  velocity 


(11-33) 
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vhere  the  state  at  any  time  is  determined  by  numerically  in¬ 
tegrating  Eqs  (11-30),  (11-31),  and  (11-32)  from  some  initial 
state. 


BRV  Equations  of  Motion 

Once  the  RV  reenters  the  earth’s  atmosphere,  somewhere 
below  120  -km  altitude,  it  begins  to  experience  an  additional 
force  due  to  aerodynamic  drag.  Even  a  MaRV  reenters  along  a 
ballistic  trajectory,  so  that  for  a  portion  of  its  flight  the 
equations  of  motion  of  a  MaRV  are  the  same  as  those  of  a  BRV. 
The  aerodynamic  drag  force  acts  opposite  to  the  air  relative 
velocity  vector  and  is  calculated  as  the  product  of  three 
quantities! 

(i)  the  dynamic  pressure,  q 

(ii)  the  nondimensional  drag  force  coefficient,  Cj^  and 

(iii)  a  reference  area,  A 

The  dynamic  pressure,  q  depends  on  the  air  density,  p  and 
the  relative  velocity  between  the  RV  and  surrounding  air, 


V-  in  the  following  manner 


(11-34) 


Choosing  A  as  the  cross  sectional  area  of  the  RV  and  defining 
the  ballistic  coefficient,  @  as 


A  /m 


(11-35) 


l!  this  analysis,  winds  are  assumed  to  be  negligible, 
so  that  air  relative  velocity  is  velocity  relative  to  the 
rotating  earth. 
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where 


m  is  the  mass  of  the  RV, 

the  acceleration  due  to  the  drag  force  is  given  by  the  vector 
equation 


(11-36) 


where 

4  is  the  unit  vector  in  the  direction  of  V  . 

V  3 

In  the  ECI  reference  frame,  the  components  of  the  drag  ac¬ 
celeration  are  determined  to  be 


(11-37) 


(11-38) 


(11-39) 


• 


where  X  ,  Y^,  and  Z  are  the  components  of  the  air  relative 

3  o  3 

velocity  in  the  ECI  reference  frame  given  by 


(11-40) 


t  --  f-m 

X 


(11-41) 

(11-42) 


•  • 


where  X,  Y,  and  Z  are  the  components  of  the  inertial  velocity 
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in  the  ECI  coordinate  system  and  X,  Y,  and  n  are  as  previously 
def ined . 

The  total  acceleration  of  the  BRV  is  the  sum  of  the 
gravity  acceleration  and  the  drag  acceleration.  So  the 
equations  of  motion  of  the  BRV  in  ECI  coordinates  are  given  as 

(11-43) 

(11-44) 

(11-45) 

•  •  •  *  •  • 

Where  X  ,  Y  ,  and  Z  are  given  by  Eqs  (11-30),  (11-31),  and 
(11-32). 

Although  BRV  flight  is  somewhat  idealized  in  that  there 
is  almost  always  some  lift  due  to  non-zero  angle  of  attack  or 
asymmetric  ablation,  the  net  effect  is  usually  assumed  to  be 
negligible  in  the  motion  model  of  high  performance  ballistic 
RV's.  The  state  vector  for  the  reentry  of  a  ballistic  RV  is 

(11-46) 

If  the  unintentional  lift  mentioned  above  becomes  significant, 
then  the  ballistic  RV  model  must  be  modified  and  will  be  the 
same  as  that  of  a  MaRV . 


laRV  Equations  of  Motion 


When  the  reentry  vehicle  undertakes  a  maneuver,  a 
third  force,  aerodynamic  lift  is  introduced.  The  lift  force 


23 


is  in  a  plane  perpendicular  to  the  air  relative  velocity. 

One  method  which  allows  the  lift  acceleration  to  be  formulated 
in  a  manner  similar  to  the  drag  acceleration  is  to  decompose 
the  lift  force  into  two  orthogonal  components.  One  component 
is  chosen  to  be  in  a  direction  whose  unit  vector,  is  given 
as  (see  Fig  II-4) 


(11-47) 


The  other  component  is  chosen  to  be  in  a  direction  whose  unit 
vector,  4^  is  given  as 


(11-48) 


These  two  components  of  the  lift  force  are  called  turn,  and 
climb,  F^.  Figure  I I -4  shows  the  geometry  of  the  lift  compo¬ 
nents  in  relationship  to  ^  and  the  drag  force,  Fj^.  With 
these  definitions,  turn  acceleration  and  climb  acceleration 
are  given  as 


(11-49) 

(11-50) 


where 

it 

is  the  non-dimensional  turn  force  coefficient,  and 

it 

is  the  non-dimensional  climb  force  coefficient 
To  simplify  the  notation,  two  new  terms  are  introduced,  , 
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the  turn  parameter,  and  C^, the  climb  parameter,  defined  below 

Cr-%:  . 

Cc  =  ^  (11-52) 

The  ECI  components  of  the  turn  and  climb  accelerations  may 
then  be  written  as 


(11-53) 

(11-54) 

(11-55) 

(11-56) 

(11-57) 

(11-58) 

(11-59) 


Now  all  the  equations  necessary  for  representing  the 
equations  of  motion  of  the  MaRV  completely  in  the  computa- 
tional-ECI  reference  frame  have  been  given.  In  some  cases 
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the  accelerations  have  been  expressed  in  terms  of  air  rela¬ 
tive  coordinates  because  of  the  convenience  of  notationi 
hovever,  in  those  cases,  expressions  relating  the  air  rela¬ 
tive  coordinates  to  inertial  coordinates  have  been  given. 
Symbolically  then,  the  equations  of  motion  for  a  MaRV  in  the 
ECI  reference  frame  are  given  as 


X  =  X.  +X,  ^Xr 

(11-60) 

#•  •• 

y=  4  *y,*yr*yr. 

(IT-61) 

^  ^  •  •  « « 

Z  -  Zj  'Zr  'Z. 

(11-62) 

•  •  •  •  •  • 

vhere  X  ,  Y  ,  and  2  are  given  by  Eqs  (11-30),  (11-31),  and 
(11-32) I  Xjj,  Yjj,  and  Zp  are  given  by  Eqs  (11-37),  (XI-38),  and 

(11-39)  I  X,j„  y^,  and  2^  are  given  by  Eqs  (II-53),  (11-54),  and 

«•  ••  •• 

(11-55) j  and  X  ,  Y  ,  and  Z_  are  given  by  Eqs  (II-56),  (11-57), 

WWW 

and  (11-58), 

Numerical  Integration  of  Equati_QDs,_of  Motion 

The  aquations  of  motion  given  by  Eqs  (11-60),  (11-61) 
and  (11-62)  represent  a  system  of  coupled,  second  order,  non¬ 
linear  ordinary  differential  equations  which  cannot  be  solved 
in  closed  form.  They  must  be  numerically  integrated  from 

•  «  e 

some  initial  conditions  on  X,  Y,  Z,  X,  Y,  and  Z  to  give  the 

/ 

position  and  velocity  at  any  other  time.  A  numerical  inte¬ 
gration  scheme  which  is  used  extcnsi^rely  in  ballistic  missile 
trajectory  simulation  and  estimation  is  a  fourth-order 
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Runge-Kutta  method  modified  by  Nystrom  to  handle  the  second 
order  equation.  Kreyszig  (Ref  21)  calls  this  the  Runge- 
Kutta-NystrSb  method.  Its  popularity  for  propagating  the 
state  of  a  ballistic  missile  is  due  to  the  fact  that  it  is 
self  starting;  it  has  small  round  off  error  for  small 
size;  the  step  size  may  be  varied  at  will  from  step  to  step; 
and  it  is  straightforward  to  implement.  The  following  de¬ 
scription  follows  essentially  that  of  (Ref  33) . 

Let  the  equations  of  motion  given  in  the  previous  sec¬ 
tion  be  represented  by  the  functional  relationship 


(11-63) 
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Furthermore,  let  r(t)  denote  the  solution  of  Eq  (11-63) 

corresponding  to  the  initial  conditions  r(t^)  =  r  and 
•  •  • 

~j  ~j  approximations  to 

and  £.(tj)  and  set  h  =  tj^^-tj.  To  obtain  the  approximations 

r  ...  and  r  ...  we  calculate  as  follows 

~j+l  -j+1 


(11-65) 
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(11-66) 

K>  --  "  a  .  rj^  i  !<•.) 

(11-67) 

(11-68) 

After  the  above  are  computed,  the  new  values 

are  then  obtained  as 

• 

and 

(11-69) 

ih  ‘  fj  z[^,  ’■  iK, 

(11-70) 

• 

Thus,  starting  with  r^  and  r^  as  Known  quantities  the  inte¬ 
gration  subroutine  which  consists  basically  of  Eqs  (11-65) 

through  (11-70),  repeats  itself  cyclically  until  as  many 

steps,  h  as  may  be  desired  have  been  completed. 


Drag  Models 


As  Athans  (Ref  3)  states,  the  ballistic  coefficient  is 
an  important  parameter  for  several  reasons.  This  parameter 
physically  represents  the  ratio  between  the  mass  of  the  RV 


and  the  effective  drag  area  along  the  velocity  vector.  Thus, 
knowledge  of  this  coefficient  reveals  some  telltale  charac¬ 
teristics  about  the  missile.  In  addition,  an  accurate  6- 
estimate  generally  imposes  a  more  stringent  requirement  than 
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do  other  criteria  on  the  estimator  performance  and,  therefore, 
offers  to  us  a  convenient  method  of  evaluating  estimator  ac¬ 
curacy  or  comparing  estimator  results. 

The  ballistic  coefficient  given  in  Eq  (11-35) 

(11-35) 

is  not  a  constant  because  the  drag  coefficient,  Cj^  varies 
with  Mach  number  for  an  RV  of  given  geometric  shape.  For 
sphere-cone  RVs,  Cj^  is  parameterized  as  functions  of  the 
vehicle  half -cone  angle,  bluntness -ratio  (radius  of  nose/ 
radius  of  base)>  and  Mach  number.  A  series  of  unpublished 
tables  from  wind  tunnel  test  of  sphere-cone  RV’s  is  stored 
internal  to  the  program  and  the  user  supplies  an  estimate  of 
the  half  cone  angle  and  bluntness  ratio.  Figure  I I -5  shows 
the  theoretical  drag  curve  reproduced  from  Ref  (35)  for  a  6® 
half -cone-angle,  0.2l-bluntness-ratio,  sphere-cone  RV. 

Since  tables  don't  lend  themselves  readily  to  analytic 
differentiation,  an  analytic  expression  for  drag  is  needed 
for  the  partials  necessary  in  the  differential  correction 
process.  The  drag  model  used  for  this  purpose  is  (Ref  34) 

where  M  is  Mach  Number  and  k^,  and  k^  are  related  to  the 
BV  half -cone  angle  and  bluntness  ratio  and  are  determined 
from  least  squares  fits  of  the  drag  tables  to  Eq  (II -71). 
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Although  Eq  (11-71)  doesn't  represent  all  Mach  Numbers  shown 

in  Fig  II-5,  it  represents  quite  adequately  the  region  in 

* 

which  ve  are  interested,  down  to  about-  Mach  3. 


Models  of  the  Atmosphere 

Models  of  the  atmospheric  density  and  temperature  are 
required  for  accurate  trajectory  simulation  and  estimation. 
The  density  is  rt;qui  red  for  calculating  the  dyneunic  pressure, 
and  temperature  is  required  for  Mach  number  calculation. 


(11-34) 


(11-73) 


The  importance  of  accurate  density  information  for  accurate 
^  calculation  is  seen  by  the  fact  that  a  percent  error  in 
density  transforms  directly-  into  a  percent  error  in  3  owing 
to  the  linear  relationship  between  P  and ^ 


iff 

m 


(11-74) 


Thus,  the  program  provides  for  the  input  of  the  local  density 
and  temperature  tables  as  functions  of  altitude.  In  addition, 
for  simulation  purposes,  tables  from  the  1962  U.S.  Standard 
Atmosphere,  are  stored  internal  to  the  program.  Again, 


*  2 

For  I CBM  ranges,  a  high  performance  RV  (P>7000Kg/m  ) 

impacts  about  Mach  3 . 
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however f  table*  don't  lend  themeelvee  to  analytic  differentia¬ 
tion.  Thue,  an  analytic  expreeeion  for  density  ie  required. 
The  density  model  which  is  used  for  the  partial*  needed  in 
the  differential  correction  process  is  (Ref  7) 


<II-75) 


where 


A 


■  (1.226Kg/m^)  the  atmospheric  density  at  h  ■  0 
e  ■  the  constant  2.7182818 


s.h.  ■  the  scale  height  7000m) 


HI  TM  gftrlfmt9r,.Bguiiifcigni 


Zn  th«  previous  chapter*  •cuatione  were  givon  for  prop¬ 
agating  tha  atata  of  a  MaRV  forward*  in  time.  Thaaa  war#  ax- 
praaaad  in  tarma  of  tha  diffarantial  aguationa  of  motion 
along  with  an  algorithm  for  numarically  integrating  tha 
aquationa  of  motion*  Additional  aquationa  ware  given  for 
tranaforming  tha  poaition  of  tha  RV  at  any  time*  t^  from  ECl 
ooordinataa*  (X^,  Y^*  Z^)  to  apharical  radar  coordinataa* 

Aa^ *  )  • 

Zn  thia  chapter*  aquationa  for  aatimating  tha  atata  of 
a  MaRV  at  aoma  epoch  time*  t^*  baaed  on  tha  obaarvad  radar 
maaauramanta  and  tha  aquationa  given  in  tha  pravioua  chapter* 
are  given*  Thaaa  aatimator  aquationa  nra  baaed  on  tha  oiaa- 
•  ioal  method  of  weighted  laaat-aquaraa*  which  ia  generally 
agreed  to  have  bean  invented  firat  by  Oauaa  for  hia  aatronom- 
ioal  atudiaa  in  179S*  Ha  waa  18  ya^ra  old  at  tha  time* 

Thera  are  two  othara*  lagandra  in  franca  in  1008  and  Adrain 
of  Amirigi  in  1M08*  Who  independently  developed  the  method. 
Since  theaa  acirly  timaa*  there  haa  bean  a  vaat  literature  on 
varioua  aapaota  of  tha  laaat-aquaraa  method  juat  as  Oauaa  pra- 
diclad  there  would  be*  for  data! la  of  tha  vary  Intareating 


Tha  aquationa  given  will  allow  for  tha  propagation  to 
be  forward  or  backward  in  time  with  equal  aaaai  however*  all 
propagationa  of  tha  state  in  this  report  will  be  forward  in 
time* 


history  of  estimation  theory,  the  reader  is  referred  to  the 
excellent  survey  papers  listed  in  Refs  (19)  and  (30)  in  the 
bibliography. 

The  basic  weighted-least-squares  differential  correc¬ 
tion  procedure  is  developed  first.  Then  the  relationship  be¬ 
tween  WLS  estimation  and  minimum  variance  and  maximum  likeli¬ 
hood  estimation  is  discvissed.  The  estimator  with  and  without 
apr.iori  statistics  is  described.  Sequential  estimation  and 
an  algorithm  for  deciding  when  to  switch  to  a  higher  order 
estimator  are  discussed.  A  method  for  obtaining  initial  con¬ 
ditions  is  also  suggested. 


lat-Souares  Differential  Correction 


The  method  of  least  squares  is  concerned  with  estima¬ 
ting  the  values  of  a  set  of  parameters  of  the  measurement 
model,  which  relates  the  measurements  to  the  parameters,  by 
minimizing  the  sum  of  the  squares  of  the  residuals.  The 
residuals  are  defined  to  be  the  differences  between  the  ob¬ 


served  measurements  and  the  corresponding  measurements  com¬ 
puted  from  the  mathematical  model  using  some  set  of  estimates 
of  the  parameters.  The  set  of  estimates  which  renders  the 
sxun  of  squared  residuals  a  minimum  is  said  to  be  optimal  in 
the  least-squares  sense.  To  put  this  in  terms  of  the  problem 
at  hand,  the  observed  measurements  are  the  ranges,  azimuths, 
and  elevations  from  the  tracking  data  and  the  corresponding 
computed  measurements  are  determined  from  the  equations  in 
Chapter  II.  That  is  to  say  that  the  equations  of  Chapter  II 
constitute  our  measurement  model.  The  parameters  we  want  to 
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estimate  are  the  RV  state  variables  at  some  fixed  time,  t^ 
called  the  epoch  time.  If  we  let  the  vector  represent 
the  set  of  observed  measurements  of  range,  azimuth,  and 
elevation  at  any  given  time,  t^,  (t^  >  t^)  the  measurement 
model  may  be  represented  notationally  as 


(iii-i) 


where 

n  =  the  total  number  of  measurement  times 

=  the  true  (but  unknown)  RV  state  vector  at  epoch 
h^  =  the  vector  of  nonlinear  functions  relating  ^  to 
=  the  vector  of  errors  in  the  observed  measurements 
at  t. 

In  component  form,  the  measurement  model  is  given  as 


(III-2) 


(III-3) 


(III-4) 


where  and  are  the  measurement  errors  in  range, 

azimuth,  and  elevation  respectively  at  t^  and  the  o  sub¬ 
script  on  ^oi'  ^oi  identify  them  as  the  ob¬ 

served  values  at  t^^ , 
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The  most  general  state  vector  for  a  MaRV  contains  nine 
elements -three  positions  elements,  three  velocity  elements, 
and  three  aerodynamic  parameters.  Thus,  the  state  vector  at 
epoch  may  be  given  as 


)(.  - (X., y., z.,  fi., Cr, Cc) 

..his  would  represent  the  MaRV  during  maneuvering  pro\'ided 
and  are  not  zero.  Ballistic-flight  is  modeled  by  setting 
and  to  zero  and  free-flight  is  modeled  by  setting  '-/Pq* 
C,j,,  and  to  zero. 

The  set  of  equations  given  by  Eqs  (III-2)  -  (III-4) 
are  the  observation  equations.  If  they  were  linear,  we  could 
taKe  any  i  of  them  and  solve  for  the  unknovms,  where  Ji  is 
the  number  of  unknown  state  variables  (i.e.,  i  *  6,  7,  or  9). 
(Here  we  have  assumed  that  i  i  3  n  )  The  problem  is  the  non¬ 
linearity  of  the  observation  equations.  There  is  no  direct 
method  for  simultaneously  solving  a  set  of  nonlinear  equa¬ 
tions;  therefore,  an  iterative  approach  must  be  taken.  The 
iterative  approach  when  applied  to  the  least-squares  observa¬ 
tion  equations  is  typically  called  differential  correction. 

It  consists  of  "linearizing"  Eq  (lll-l)  by  expanding  it  in  a 
Taylor's  series  about  some  initial  guess  or  estimate  of  X^j 
truncating  it  by  ignoring  all  terms  of  order  two  or  higher; 
and  solving  for  the  resulting  differentials  which  minimize 
the  sum  of  the  squared  residuals  determined  from  using  the 
initial  estimate.  This  optimal  set  or  a  fferentials  is  then 
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0^, 


added  to  the  initial  estimate  to  form,-  a  better  estimate,  and 

the  process  is  repeated  until  some  convergence  criterion  is 

met  (e.g.,  all  of  the  differentials  become  essentially  zero). 

This  is  explained  in  more  detail  in  vhat  follows  where  the 

weighted -least-squares  estimate  is  derived. 

A 

Let  X,  be  an  estimate  of  3^,  then  the  Taylor’s  series 
— o  .  — o  ■' 

expansion  of  Eqs  (III-2)  -  (III-4)  are  given  as 


(III-6) 


(III-7) 


(III-8) 


where  the  vector  of  differentials,  ^  for  the  free-flight 
case  is  given  as 


V.-9. 

X-L 


(III-9) 


The  partials  of  range,  azimuth,  and  elevation  with  respect  to 
the  state  vector  ^  given  in  Eqs  (III -6)  -  (III-8)  are 
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A 

evaluated  using  the  current  estimate,  ^  and  are  defined  by 


the  following  matrix  equation 

IS  iS  U  'i£ 

3A.  ay.  az.  ^  3%  az, 

^  ^  t  ^  ^ 

3.  ^  tj, 

The  set  of  equations  given  by  Eqs  (III -6)  -  (III -8)  when  the 
higher  order  terms  are  neglected  is  a  set  of  linear  equations 
in  the  unknowns  ^X*  is  this  set  of  equations  that  we 

apply  the  principle  of  least-squares.  First  we  make  the 
notationally  simplifying  definitions 


(III-ll) 


(III-12) 


Some  of  the  reelduale  of  Bq  (111-14)  are  in  terms  of  range 
^  units  and  some  are  in  terms  of  angle  units i  and  even  though 

the  range  measurements  are  more  accurate  than  the  angle  meas¬ 
urements,  the  numerical  values  of  the  range  residuals  (in 
meters)  will  most  certainly  be  larger  than  the  angle  residuals 
(in  degrees).  This  forces  us  to  apply  weights  to  the  observa¬ 
tions  to  normali;^Q  the  residuals  so  that  all  are  given  equal 
consideration.  Or  ve  might  want  to  apply  weights  even  if  all 
the  observations  were  of  one  type  (say  range)  to  give  more 
accurate  observations  greater  influence  and  leas  accurate 
observations  less  influence.  Assuming  we  may  want  to  assi^ 
a  distinct  weight  to  each  observation  ve  define 

Wpi  "  weight  applied  to  the  i^^  range  observation 
"  weight  applied  to  the  i^^  azimuth  observation 
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Wgi  *  weight  applied  to  the  ‘  elevation  observation 

and 

*=  the  3n-by-3n  diagonal  weighting  matrix  whose 

diagonal  elements  are  '^Ei*  ^  “  1*2,.., n. 

With  the  above  definitions,  the  weighted  observation  equa¬ 
tions  are  given  by  premultiplying  Eq  (III-14)  by  W**. 


(III-15) 


Where  the  subscript  W  on  ^  is  to  indicate  that  this  is  the 
weighted  error  function. 

Applying  the  principle  of  least  squares  to  Eq  (III-15)  we 
obtain  the  sum  of  the  squared  residuals  to  be  - 


(III-16) 


Taking  the  partials  of  the  above  scalar  function  with  respect 
to  4X  gives 


(III-17) 


Setting  Eq  (III-17)  to  zero,  transposing  and  solving  for  ^ 
results  in  the  standard  weighted-least-squares  estimate 


•  (h"'wh)'h'wV 


(III-18) 
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V/* 


(III-19) 


where 


Had  the  Eq  (III-.2)  -  (III-4)  been  linear,  the  optimal  least- 

squares  estimate  of  would  be  given  by  adding  the  ^  de- 

A 

termined  from  Eq  (III-18)  to  the  initial  estimate,  5^. 

As  it  is,  however,  since  the  observation  equations  are  not 

linear,  the  optimal  estimate  must  be  determined  in  stages  by 

adding  the  differential,  ^  determined  from  Eq  (III-18)  at 

A 

each  stage  to  the  current  estimate,  3^  and  iterating  the  above 
equation  until  the  process  converges  to  the  optimal  weighted- 
least-squares  estimate. 

In  algorithm  form  the  weighted-least-squares-dif ferential- 

correction  method  may  be  summarized  by  the  following  steps i 

A 

1.  Somehow  determine  initial  estimates,  X  for  all  of 

the  unknown  parameters  (e.g.  see  page  54) . 

A  A 

2.  Using  the  estimate,  3^,  propogate  ^  in  time 
(integrate  to  the  time  of  each  observation  and 
evaluate  V  and  H, 

3.  Determine  ^  from  Eq  (III-18). 

A 

4.  Add  AX  to  X  to  form  a  better  estimate 

—  — o 


A 

X, 


X  +  AX 
— o  — 


5.  Return  to  Step  2,  and  repeat  Setps  2-4  checking  for 
★ 

convergence  after  Step  4  of  each  iteration. 


As  given  on  page  38, 
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6.  At  convergence  the  optimal  least-squares  estimate 
A 

is  (Notei  no  notation  distinction  is  made  be¬ 

tween  the  converged  estimate  and  the  initial  est- 
A 

mate,  3^) 

In  step  4  the<=i denotes  replacement  or  writing  over. 

The  above  procedure  converges  quite  rapidly  (<  4  itera¬ 
tions)  for  good  initial  estimates.  With  poor  initial  esti¬ 
mates,  the  above  procedure  may  not  converge  at  all.  In  fact, 
it  will  most  likely  diverge.  A  technique  for  determining 
good  initial  estimates  is  given  later  in  this  chapter. 

Minimum  Variance  and  Maximum  Likelihood  Estimation 

The  weighted-least-squares  estimator  developed  in  the 
previous  section  was  done  so  from  purely  "reasonableness" 
arguments.  No  knowledge  of  the  underlying  probability  dis¬ 
tributions  of  the  observation  errors  6^.)  was  re¬ 

quired  and  none  was  used.  Weights  were  assigned  in  a  some¬ 
what  arbitrary  manner  where  it  seems  reasonable  that  more 
accurate  measurements  should  receive  nore  weight  than  less 
accurate  measurements.  As  a  result  nothing  can  be  said  about 
the  accuracy  of  the  WLS  estimate.  Clearly  from  Eq  (III-18), 
the  least  squares  estimate  is  a  function  of  W.  Therefore,  it 
appears  that  there  may  be  an  ''optimal"  W  to  use  in  the  least- 
squares  estimator.  Junkins  (Ref  18)  shows  that  if  the  opti- 

A 

mality  criterion  is  to  determine  the  estimate  ^  with  mini¬ 
mum  covariance  matrix,  then  W  should  be  the  inverse  of  the 
covariance  matrix  of  the  observation  errors.  When  the 
weighting  matrix  is  so  defined,  the  least-squares  estimator. 
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in  the  linear  case,  is  called  the  minimum  variance  (MV) 

estimator.  Defining  as  the  covariance  of  the  observation 

A 

errors  and  P  e  the  covariance  of  3^,  the  first  order  approx 
imation  to  a  MV  estimator  for  our  nonlinear  problem  is 
given  as 


41.  •  (.H\" H)'' (III-20) 

P  •  (111-21) 

A 

Where  3^  is  determined  exactly  as  in  Steps  1-6  of  the  previous 

section  with  W  now  replaced  byQ”^  and  P  is  determined  using 

i\ 

the  converged  Since  P  of  Eq  (III-?  )  is  only  a  first- 

order  approximation  of  the  covariance  matrix,  it  may,  for 

very  nonlinear  problems,  be  a  poor  estimate  of  the  errors  in 
A 

iC,  Therefore,  Monte  Carlo  simulations  are  usually  required 
to  validate  the  covariance  matrix.  This  is  discussed  more 
fully  in  Chapter  IV. 

Since  the  diagonal  elements  of  the  covariance  matrix, 

P  are  the  mean  squared  errors  of  the  individual  state  vari¬ 
ables,  the  MV  estimator  given  by  Eqs  (III-20)  and  (III-21)  is 
sometimes  referred  to  as  the  minimum  mean  square  (MMSE)  esti¬ 
mator.  See  Liebelt  (Ref  23).  It  is  also  (Ref  32)  called 
the  Markov  estimator . 

Another  assumption,  which  is  almost  always  made  in 
actual  application,  is  that  the  observation  errors  are 
Gaussian  distributed  with  zero  mean.  Using  this  assumption, 
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Wiesel  (Ref  36)  develops  the  maximum  likelihood  (ML)  esti¬ 
mator  which  is  identical  to  Eqs  (III-20)  and  (III-21)  where 
Q  is  now  the  covariance  of  the  normal  distribution  of  the 
-observation  errors.  For  ML  estimation,  see  also  (Ref  16,18,24). 

Another  assumption  which  was  made  in  developing  the 
estimator  for  this  application  is  that  the  observational 
errors  are  not  correlate.^  in  time  nor  spatially.  This  is 
not  strictly  true  spatially  since  the  radar  doesn’t  make 
its  measurements  in  range,  azimuth,  and  elevation  but  in¬ 
stead  transforms  "raw”  measurements  into  range,  azimuth,  and 
elevation.  In  addition,  if  the  measurements  are  from  a 
moving  sensor  which  has  been  transformed  to  a  fixed  location, 
then  the  errors  are  correlated  in  time.  The  effect  of  this 
approximation  is  to  make <5  a  diagonal  matrix  with  the  error 
variances  on  the  diagonal.  Considerable  storage  and  calcula¬ 
tions  can  be  saved  if  the  error  variances  are  identical 
from  time  to  time,  i.e.. 


< 

1) 

(III-22) 

< 

II 

(III-23) 

< 

(III-24) 

.  2  2  2 
for  any  two  times,  t.  and  t.  where  and  a_.  are  the 

J.  J  Kl  AX  t,X 

variances  in  range,  azimuth  and  elevation  at  t^.  This  being 
the  case,  the  3n-by-3n  covariance  matrix,  can  be  represented 
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by  three  scalars,  and  Og.  The  a^,  a^,  and  used 

in  the  numerical  examples  to  validate  the  estimator  equations 
and  the  programming  are  given  in  Table  I-l  of  Chapter  I. 

A  Priori  Statistics 

It  was  mentioned  earlier  that  good  initial  estimates 
of  the  unknown  state  parameters  are  required  for  the  differ¬ 
ential  correction  process  to  converge.  If  in  addition, 
statistics  for  the  initial  estimates  are  also  known  the  esti¬ 
mator  is  modified  to  reflect  these  a  priori  statist.’  cs  (Ref  2) 
as 

42  *  (III-25) 

P  «  (III -26) 

A 

Where  P_  is  the  covariance  of  the  initial  estimate, 

o  '  ^o 

The  Sequential  Estimator 

The  weighted-least-squares  estimator  developed  in  the 
preceding  sections  assumed  a  deterministic  system  model 
where  the  system  model  is  given  as 

X  =/Cx,y,  /,  X,  y;  Cr,Cc ,  t)  (111-27) 

Another  way  of  saying  this  is,  if  the  unknown  state  variables 
were  known,  we  could  propagate  the  state  forward  in  time 
perfectly.  That  would  imply  the  models  for  gravity,  drag, 


and  lift  given  in  Chapter  ZI  are  exact.  Hhereaa  no  model 
ie  ever  perfect#  some  are  very  nearly  so  such  that  determin¬ 
istic  least  squares  vorKs  exceptionally  veil.  A  case  in 
point  is  the  missile  in  free-flight  vhere  the  only  force 
acting  is  gravity.  Gravity  can  be  modeled  vith  sufficient 
accuracy#  and  is  in  this  algorithm#  that  all  of  th#  free- 
flight  data  can  be  processed  simultanaously  in  a  single  batch 
using  the  six-state  estimator  given  previously.  On  the  other 
hand#  the  drag  model  for  an  RV  is  not  so  veil  knovn#  and  the 
lift  model  is  even  leas  veil  knovn.  This  motivates  us  to 
seek  an  estimator  for  the  reentry  flight  regime  to  process 
the  data  sequentially#  a  small  batch  at  a  time.  The  sise  of 
tho  batch#  specified  by  the  number  of  time  points#  should#  as 
a  general  rule#  be  the  largest  over  vhich  the  model  may  be 
considered  perfect  (sufficiently  accurate).  This  number 
varies  vith  altitude#  accuracy  of  the  observation  data#  etc. 
and  no  attempt  vas  made  in  this  research  to  determine  the 
optimum  batch  siae.  Instead  the  user  of  the  program  is  al- 
lovsd  to  specify  the  batch  sise  anyvhere  from  one  point ‘to 
all  the  data.  Examples  win  be  given  in  Chapter  ZV. 

The  sequential  estimator  vhioh  was  developed  for  this 
problem  is  based  on  )7qs  (III-25)  and  (XZI-26)  and  an  addi¬ 
tional  equation  for  propagating  the  covariance  in  time  (Ref  19). 

( n  I  -2 0 ) 
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Th«  ••qiMntikl  Mtimator  vhicAi  vaa  Implatnontad  la  aimilar  co 
that  auggaatad  by  Barkar  (Raf  4}  and  ia  outlinad  in  tha  fol- 
loving  atapai 

1*  Tha  numbar  of  peinta  in  batch  ona  (NPBl*)  ia 

A 

apaoifiadi  and  tha  initial  aatimataCf  and 
ara  dotarninad*  (Tha  aubaeript,  1  rafara  to  tha 
toateh  numbar). 

2.  Tha  MLS  diffarantial  oorraution  algorithm  (atapa 
2-5)  givan  pravioualy  ia  itaratad  until  convarganca 
la  aehiavad  uaing  Bg  (ZXX-2S)  inataad  of  (IZX-IS)# 
i  «a .  f 

ii>  •  m (jij-29) 

3*  Aftar  atap  2  haa  gonvargad»  updata  tha  oovarianoa 
uaing  Bq  (ZXZ-26) 

I 

/f  ■  (xxi-30) 

'  \ 

4 1  Propagata  tha  optimal  aatimittf  data^minad  for 
A 

batch  ona 9  oovarianoa  P^  to  tha  tima  of 

tha  next  apooh.  (Oanarally  apaaking  the  apooh  ia 
advanoad  ona  tima  intar val  but  thia  ia  aalaotad  by 
tha  uaar  by  apaoifying  tha  numbar  of  pointa 


* 


Program  input 


between  epochs  one  and  two,  NPBE12*) 

(III-31) 

^  (III-32) 

(The  number  of  points  in  batch  two  and  succeeding 
batches  is  specified  separately  as  NPBAl  so  that 
batch  one  may  Incorporate  all  of  the  free  flight 
data  if  desired) 

5.  Again  the  WLS  differential  correction  algorithm  is 
iterated  until  the  optimal  estimate  of  ^2  con¬ 
verged  upon* 

A 

6«  After  stop  5  has  converged,  the  covariance  of  5^02 
is  updated  by 

Fi  •  (,Si -thf^ny  (111-33) 

A 

7#  Propagate  the  optimal  state  estimate,  ^2 

variance  P2  to  the  time  of  the  next  epoch  where 
the  time  of  the  next  epoch  is  specified, by  the 
number  of  points  between  epochs,  NPBE*.  (NPBE  and 
NPBE12  are  specified  separately  again  so  that  batch 
one  may  be  treated  separately). 


*Program  input,  NPBE12  <  NPBl  and  NPBE  <  NPBAl 
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6. 


Repeat  steps  5  and  6  for  batch  3. 

9.  This  propogate-update  cycle  is  repeated  for  sub¬ 
sequent  batches  until  the  final  data  point  is  in¬ 
cluded  in  a  batch. 

10.  Once  the  final  point  is  included  in  a  batch,  the 

batch  size  is  gradually  collapsed  to  a  single  data 
point . 

Thus,  the  sequential,  weighted-least-squares  estimator  is 
controlled  by  four  inputs  NPBl,  NPBE12,  NPBAl  and  NPBE 
where 

NPBl  =  number  of  points  in  batch  one 
NPB12  =  number  of  time  intervals  between  epochs  one 
and  two 

NPBAl  =  number  of  points  in  all  batches  after  one 
NPBE  =  number  of  time  intervals  between  succeeding 
batches 

When  NPBl  is  set  to  3n  or  greater,  all  of  the  data  are  pro¬ 
cessed  simultaneously  in  a  single  batch.  When  NPBl  and  NPBAl 
are  both  set  to  one  and  the  differential  correction  algorithm 
is  not  iterated  (MAXIT  =  0),  the  equations  make  up  one  form 
of  the  Kalman  filter  (see  Ref  24). 

Finite  Memory  Estimator 

Experience  gained  in  developing  the  estimator,  indicates 
that  for  data  simulated  as  indicated  in  Table  I  if  the  batch 
sizes  are  50  points  or  greater,  then  no  a  priori  statistics 
are  needed  (i.e.,  P”^  =  0,  all  i).  On  the  other  hand,  when 


the  batch  size  is  sufficiently  small  that  the  normal  matrix, 

(H  Q  H)  is  ill-conditioned  or  singular,  a  priori  statistics 
are  required.  In  this  small  batch  size  applications,  it 
has  been  found  that  with  the  infinite  memory  estimator  given 
above  which  includes  Eq  (III-26),  that  the  covariance  matrix 

P  '  ip"* (III-33) 

becomes  indefinite  as  the  batch  number,  i  becomes  larger.  To 

★ 

prevent  this  from  happening,  a  technique  described  by  others 
was  implemented  to  deweight  the  a  priori  covariance.  This 
results  in  modifying  Eq  III-33  as 

p  ’  ic p" * H'a' H)"  (III-34) 

where  C  is  a  diagonal  matrix  (Barker  used  a  scalar)  whose 
diagonal  elements  are  between  zero  and  one.  Many  simula¬ 
tions  were  run  to  determine  a  C  to  keep  positive  definite. 
No  conclusions  can  be  drawn  but  it  seems  that  the  larger  the 
batch  size,  the  smaller  C  must  be.  Small  batch  sizes  re¬ 
quired  >  .8.  One  way  to  avoid  this  problem  would  be  to 
use  large  batch  size  (NPBAl  >50)  and  set  =  0  and  this 
appears  to  work  fine  for  ballistic  flight;  however,  if  the 
MaRV  undergoes  some  high-g  maneuvers,  small  batch  sizes  would 
be  necessary  as  the  coefficients  of  turn,  and  climb, 
are  modeled  as  constant  over  the  batch  time  interval .  The 
case  of  high-g  maneuvering  was  not  investigated  in  this  study. 

*This  deweighting  accomplished  in  various  ways.  See 
(Ref  4,16,14).  51 


Adaptively  Increasing  the  Estimator  Order 

As  has  been  mentioned  previously,  during  free-flight, 
the  estimator  is  a  six-state  estimator  corresponding  to  the 

•  •  •  Ip 

six  elements  of  the  free-flight  state  vector,  (X,Y,Z,X,Y,Z)  . 
In  the  sequential  mode,  as  the  estimator  is  stepping  through 
the  data,  there  is  a  point  after  which  the  free-flight  esti¬ 
mator  is  no  longer  adequate.  This  is  a  result  of  the  slow¬ 
down  of  the  RV  due  to  drag  making  the  free-flight  equations 
no  longer  adequate  for  modeling  the  RV’s  motion.  At  this 
time,  we  want  to  switch  to  the  seven-state  BRV  estimator  de- 

...  f 

fined  by  the  BRV  state  vector  (X,Y,Z,X, Y,Z,  .  However, 
this  point  is  not  well  defined  a  priori.  If  we  switch  to 
the  seven-state  estimator  too  soon  (i.e.,  before  the  ballis¬ 
tic  coefficient  becomes  observable),  then  the  normal  matrix, 
(H  Q  H)  will  be  singular,  and  therefore  not  invertible. 

This  condition  will  lead  to  the  estimator  failing.  The 
penalty  for  switching  late  is  not  so  great.  It  will  just 
mean  that  for  a  period  of  time,  the  RV  motion  is  modeled  in¬ 
accurately.  This  inaccurate  modeling  of  the  RV  motion  will 
lead  to  larger  than  expected  residuals.  This  divergence  of 
the  residuals  is  what  is  used  to  determine  when  to  switch  to 
the  higher  order  estimator.  This  is  accomplished  by  keeping 

track  of  what  is  called  the  total  root  mean  square,  (RMST) 

* 

error  of  the  range  residuals 


ir 

Originally  azimuth  and  elevation  RMS's  were  also  used 
but  too  many  false  switches  occurred.  See  Example  2  Chapter  IV. 
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RMsr  =1 


Al-I 


(III-35) 


where 

=  the  range  residual 

N  =  the  total  nvunber  of  range  observations  which 
have  been  processed  at  any  given  time 
For  N  greater  than  30,  the  RMST  should  be  very  nearly  equal 
to  the  standard  deviation  of  the  normal  distribution  of  the 
range  errors,  which  for  the  simulated  data  is  five  meters. 

A  long  free-f light  track  is  useful  in  determining  a  good 
value  of  the  RMST. 

The  RMST  is  tised  to  detect  divergence  by  comparing  it 
to  the  range  RMS  for  any  batch.  Call  the  range  RMS  calcu¬ 
lated  from  a  single  given  batch  the  BRMS.  In  the  limit, 
BRMSj^-RMSTj^  Should  be  less  than  2RSTDj^  with  a  probability  of 
about  .95,  so  if 

SL/^ST7>^  (IIJ-3C) 

then  divergence  has  occurred  and  the  switch  is  made,  where 
RSTD  is  a  measure  of  the  dispersion  of  the  BRMS's  defined  as 

(III-37) 


where  K  is  the  number  of  batches  processed  thus  far.  The 
check  is  not  made  until  K  is  greater  than  four  to  allow  time 
for  "settling"  of  RSTD.  The  value  four  is  arbitrary. 
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Onc^  the  svitch  from  the  six -state  estimator  to  the 
seven-state  estimator  has  loeen  made,  the  last  batch  is  re¬ 
processed  with  the  higher  order  estimator*  and  then  the  same 
logic  is  used  to  detect  the  onset  of  maneuvering*  If  a  man¬ 
euver  is  detected*  the  program  switches  to  the  nine-state 

•  •  •  ff 

estimator  given  by  the  MaRV  state  vector*  (X*Y*?.*X*Y*Z,0 *c^c^) 

Oood  initial  conditions  for  the  state  variables  are 
required  for  the  differential  correction  process  to  converge. 
The  state  variables*  ^*  C^,  and  are  not  so  critical  and  a 
rough  guess  ie  sufficient*  e.g**  ^  •  70po  kg/m  *  ■  0  and 

■  0.  On  tho  other  hand*  the  position  and  velocity 
initial  conditions  are  very  critical  and  there  is  no  way  to 
guess  them*  Thus  a  routine  vas  written  to  determine  initial 
estimates  as  follows 

1.  The  range*  azimuth*  and  elevation  measurements  are 
transformed  into  BCI  X*  Y*  and  Z, 

2.  The  BCI  coordinates  are  fit  in  a  lesst-squartes  sense 
to  a  special  4  order  polynomial  in  time  (given 
here  only  the  X  coordinate) 

Xj  *a,  ■* t  (III-36) 

i  • 

where  a^*  a^  and  a^  are  first  determined  from 

6 dfitrt,)  )  I*  >,%,••• 


/ 


(III-39) 


where  is  determined  froin  the  gravity  subroutine 
using  the  raw  X,Y,2  coordinates. 

3.  Once  the  a's  have  been  determined,  the  fourth  order 

polynomial  is  differentiated  to  determine  velocity. 

« 

If  all  the  data  are  free-f light,  the  above  procedure 
works  exceptionally  well  resulting  in  convergence  of  the  dif¬ 
ferential  correction  algorithm  in  one  or  two  iterations.  It 
also  works  very  well  with  data  which  includes  reentry  pro¬ 
vided  a  substantial  portion  is  free  flight.  If  all  the  data 
are  reentry,  better  initial  conditions  are  provided  by  ig¬ 
noring  the  second-order  polynomial  constraint  on  a^,  and 
Sg  given  in  step  2  i.e.,  Eq  (III -39). 

* 

This  technique  for  determining  initial  position  and 
velocity  was  developed  by  the  author  originally  for  program 
TEAP,  a  program  used  at  FTD  for  fitting  free-f light  data. 
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IV  Numerical  Exa 


moles 

Many  numerical  simulations  were  run  to  validate  the 
estimator  equations  and  the  programming  of  them.  The  pro¬ 
gram  was  designed  in  stages  with  each  stage  being  checked  for 
accuracy  before  proceeding  to  the  next  stage.  The  first 
stage  consisted  of  the  free-flight  equations  of  motion  with 
the  batch  WLS  estimator  equations.  Numerical  partials  were 
generated  to  check  the  accuracy  of  the  analytic  partials 
used  in  the  estimator  equations .  Subsequent  stages  included 
the  BRV  equations  of  motion,  the  sequential  estimator  equa¬ 
tions,  and  finally  the  MaRV  equations  of  motion  with  the  re¬ 
quired  partials.  Numerical  partials  were  simulated  for  each 
stage  to  test  the  accuracy  of  the  analytic  partials .  These 
tests  verified  the  fact  that  the  simplified  equations  of 
motions  used  for  the  partials  calculations  were  of  suffi¬ 
cient  accuracy. 

Once  the  program  was  completed,  many  more  simulations 
were  run,  and  refinements  and  modifications  were  made  where 
necessary.  The  results  of  several  of  these  have  been  al¬ 
luded  to  in  previous  chapters.  It  was  found,  for  instance, 
that  the  a  priori  covariance  matrix  had  to  be  deweighted  to 
maintain  positive  variances  as  output  of  the  estimator  equa¬ 
tions.  Simulations  showed  that  the  deweighting  had  to  be 
increased  as  the  batch  size  increased.  Because  of  the  time 
involved  in  tuning  the  estimator  equations  for  various  batch 


56 


sizes,  no  conclusion  can  be  reached  as  to  how  best  to  use 
the  a  priori  covariance  information.  It  was  determined, 

however,  that  for  batch  sizes  of  10  or  more  points  no 

.  .  .  T  -1  .  . 

a  prion  covariance  was  necessary  to  keep  HQ  H  positive 

definite.  As  a  result,  the  numerical  simulations  discussed 

in  the  remainder  of  this  chapter  were  run  with  no  a  priori 

covariance. 

Two  numerical  examples  are  presented  in  detail  in  the 
remainder  of  this  chapter.  The  first  discusses  the  opera¬ 
tion  of  the  estimator  in  the  batch  WLS  mode  where  all  data 
are  fit  simultaneously.  The  second  example  discusses  the 
results  of  running  the  estimator  in  the  sequential  mode.  A 
ten  run  Monte  Carlu  simulation  was  performed  for  each  of  the 
two  examples . 

Example  One  -  The  Batch  Estimator 

The  purpose  of  this  example  was  twofold,  (1)  to  verify 
the  accuracy  of  the  models  used  in  the  program  and  (2)  to 
look  at  the  accuracy  that  could  be  attained  when  using  the 
estimator  in  the  batch  WLS  mode  where  all  the  data  are  fit 
simultaneously.  This  mode  represents  the  ultimate  in  the 
accuracy  which  can  be  attained  from  radar  data  of  a  given 
accuracy  and  is  applicable  to  the  RV  in  free-flight  and  to 
the  RV  for  ballistic  reentry.  In  order  to  accomplish  the 
first  goal,  the  data  used  in  this  example  were  simulated 
using  the  Foreign  Technology  Division's  MRVRAP  program  -  a 
program  which  has  been  tested  and  proven  t.  >  be  accurat  e  for 
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Bimulating  the  flight  of  a  ballistic  RV.  Ten  sets  of  simu¬ 
lated  radar  tracking  data  were  generated.  Each  set  of  data 
was  corrupted  by  random  noise  from  normal  distributions  having 
zero  means  and  standard  deviations  of  5  m,  0.015  deg,  and 
0.020  deg  for  range,  azimuth,  and  elevation,  respectively. 

The  errors  in  the  data  were  determined  from  a  random  number 
generator  by  specifying  the  ••SEED*'t  The  SEED  was  specified 
to  be  different  for  each  data  set,  so  that,  even  though  the 
errors  for  each  set  were  from  the  same  three  normal  distribu¬ 
tions,  the  samples  wert  different  for  each. 

Data  for  this  first  example  were  simulated  for  a 
ballistic  RV  having  a  hypersonic  ballistic  coefficient  of 
3300  kg/sq  m.  The  ballistic  coefficient  was  modeled  to  vary 
with  Mach  number  as  discussed  in  Chapter  II.  Figure  IV-1  is 
a  plot  of  normalized  ballistic  coefficient,  drag  acceleration, 
air  relative  velocity  and  altitude  plotted  versus  time  from 
120  km  to  impact  to  show  how  these  parameters  are  related  to 
one  another  through  out  reentry.  Epoch  was  defined  to  be  at 
reentry  (120  km  altitude)  and  only  the  data  from  reentry  to 
impact  were  considered.  Figure  IV-2  is  a  plot  of  one  sample 
of  the  trackirjg  data  used.  The  parameters  have  been  normal¬ 
ized  so  as  to  be  able  to  display  them  on  a  single  graph. 

Ten  points  per  second  were  simulated  for  a  total  of  500  data 
points . 

All  of  the  data  in  each  set  were  fit  simultaneously 
solving  for  the  nine  elements  of  the  MaRV  state  vector.  The 
initial  conditions  for  the  position  and  velocity  states  were 

if 

The  SEED  is  just  a  starting  value  for  the  random 
number  generator. 
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determined  automatically  by  the  program  ae  described  in 

Chapter  III)  vas  initialised  at  500  kg/aq  m  and  both  C,j, 

and  vere  initialized  at  0*0 • 
c 

Since  the  data  for  this  example  vere  generated  from 
known  conditions,  the  performance  of  the  estimator  can  be 
judged  by  hov  veil  the  estimated  parameters  match  the  true 
state  parameters.  Table  IV-l  is  a  list  of  the  results  from 
the  ten  simulations.  In  the  upper  left  block  of  Table  lV-1 

A 

are  given  the  true  errors  in  the  state  estimates  for 

the  ten  runs.  The  last  rove  of  Table  lV-1  give  the  means 
and  standard  deviations  for  the  ten  run  sample;  The  last 
three  rows  of  the  tablo,  give  the  standard  deviations  of  the 
range  residuals,  o^.  All  runs  converged  in  three  Iterations. 
Table  IV-2  contains  the  standard  deviation-correlation  coef¬ 
ficient  matrix  from  one  of  the  single  sample  runs.  Since 
this  matrix  is  symmstric,  only  the  lovsr  triangular  form  ia 
given.  This  matrix  ia  determinad  from  the  estimator  oo- 
vsrisncs  matrix,  .  Ths  diagonal  elemsnts  are  ths 

standard  deviations.  Comparing  the  standard  dovisticns  from 
Tabls  IV-2  with  ths  last  column  of  Table  XV-1  shows  s  vary 
good  agrosmsnt . 

Anothsr  indication  as  to  hov  well  the  model  agrees 
with  the  data  is  givsn  by  plotting  the  residuals,  yigurs 
lV-3  in  a  plot  of  ths  range,  ssimuUi,  and  elevation  rosidunls 
from  one  of  the  single  sample  runs.  The  absence  of  trend  in 
the  reoiduals  Indicates  that:  the  data  have  boon  modeled 
correctly  and  accuratoly. 
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covariance  matrix 
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iwo  -  me  aeguenriai  estimator 

This  example  was  designed  to  show  how  the  sequential 
mode  of  the  estimator  operated  as  the  number  of  points  per 
batch  was  varied.  Three  specific  things  were  of  interest i 

(1)  how  the  estimator  switching  algorithm  operated  as 
the  number  of  points  per  batch  was  varied, 

(2)  how  the  accuracy  in  the  estimated  parameters 
varied  with  the  number  of  points  per  batch,  and 

(3)  how  well  the  estimator  covariance  tracked  the  true 
errors  in  the  estimated  parameters. 

For  this  example  another  ten  sets  of  noisy  data  were  gener¬ 
ated.  All  ten  samples  were  identical  except  for  the  random 
noise  added  to  the  data  which  was  varied  as  was  explained  in 
Example  One.  The  data  were  simulated  this  time  starting  at 
about  150  seconds  prior  to  reentry  and  going  to  impact  for  a 
total  of  206,5  seconds.  The  data  were  simulated  as  1  pps 
prior  to  reentry  and  10  pps  thereafter.  This  amounted  to  a 
total  of  698  time  points,  and  as  in  the  previous  example, 
the  errors  in  the  data  were  simulated  as  zero  mean,  white 
Gaussian  noises  with  standard  deviations  of  5  m,  0.015  deg, 
and  0.020  deg  respectively  for  range,  azimuth,  and  elevation. 
One  of  the  data  sets  has  been  plotted  in  Figure  IV-4  for  re¬ 
ference.  The  data  were  simulated  with  a  hypersonic  ballistic 
coefficient  of  10,000  kg/sq  m  and  5  varying  with  Mach  number. 
A  moderate  maneuver  was  begun  at  approximately  10  seconds 
prior  to  impact  (~  25  km  altitude)  by  setting  to  0.000001 
snd  to  0.00001  sq  m/kg  at  this  time.  These  values  were 
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selected  to  extend  the  flight  time  by  about  five  seconds. 

It  was  found  that  a  value  of  much  larger  than  the  one 
used,  applied  at  this  altitude,  resulted  in  the  RV  skipping 
out.  The  plus  sign  on  implies  a  left  turn  and  on 
implies  a  climb  up.  The  pertinent  reentry  parameters  of 
this  model  are  plotted  in  Figure  IV-5. 

Three  separate  cases  were  run.  In  each  case  all  of  the 
free  flight  data  (first  152  seconds)  were  included  in  batch 
one.  Then  the  epoch  was  propagated  forward  150  seconds  to 
the  time  of  reentry.  Thereafter  each  batch  consisted  of 
50  points  for  case  one,  30  points  for  case  two  and  10  points 
for  case  three  with  the  epoch  being  advanced  2  points  be¬ 
tween  batches.  Initial  estimates  of  position  and  velocity 
at  the  first  epoch  were  determined  automatically  by  the  pro¬ 
gram  and  ballistic  coefficient,  C^,  and  wore  initialized 
at  5000  kg/sq  m,  0.0,  and  0.0  sq  mAg  respectively.  The  pro¬ 
pagated  state  estimates  were  used  as  initial  estimates  for  . 
subsequent  batches.  For  each  case,  the  estimator  was  ini¬ 
tialized  as  the  six  state  free-flight  estimator  consisting 
only  of  the  position  and  velocity  states  with  the  program 
switching  to  the  BRV  and  MaRV  estimators  when  divergence  of 
the  range  residuals  was  detected.  This  shows  up  vividly 
in  Figure  IV-6  where  for  case  one,  the  total  RMS  (RMST)  of 
the  range  residuals  and  the  RMS  of  the  range  residuals  for 
each  individual  batch  (BRMS)  are  plotted  as  a  function  of 
time.  The  time  on  this  figure  is  in  seconds  after  the  first 
data  point.  Recall  that  the  estimator  switches  from  six 


states  to  seven  states  when  the  magnitude  of  the  BRMS  be- 
cones  significantly  larger,  in  a  statistical  sense,  than 
the  RMST,  This  occurs  at  approximately  182  seconds  in 
Figure  IV-6  when  the  difference,  BRMS-RMST,  becomes  larger 
than  twice  the  standard  deviation  of  the  error  in  RMST.  A 
check  is  also  made  on  the  value  of  the  dynamic  pressure  at 
this  time  to  prevent  the  filter  from  switching  while  aero¬ 
dynamic  parameters  would  still  be  unobservable.  This  pre¬ 
vented  the  filter  in  this  case  from  switching  at  3  earlier 
times  when  the  BRMS  went  above  the  RMST  +  2a  curve.  The 
estimator  switches  again  at  approximately  190  seconds  when 
the  maneuver  is  detected,  this  time  switching  to  the  nine- 
state  MaRV  estimator.  The  RMST  goes  straight  line  after 
this  switch  because  it  is  no  longer  updated  once  the  switch 
has  been  made  to  the  MaRV  estimator.  No  provision  was  made 
to  switch  from  a  higher  to  lower  order  estimator,  but  then 
it's  not  considered  necessary  since  once  the  switch  is  made 
C,p  and  will  be  observable  thereafter  even  though  they  may 

be  zero.  Also  plotted  on  Figure  IV-6  is  the  mean  of  the 
range  residuals  for  each  batch.  Ideally  this  should  be  zero 
for  each  batch.  The  same  curves  for  the  azimuth  and  eleva¬ 
tion  residuals  are  plotted  as  Figure  IV-7  and  IV-8  and  as 
can  be  seen  for  this  case  neither  azimuth  nor  elevation  de¬ 
tected  the  divergence.  This  varied  however  with  data  sets. 
In  some  the  maneuver  showed  up  dramatically  in  the  elevation 
residuals  (see  for  example  Figure  IV-9) .  In  every  case-one 
run  however,  the  divergences  were  detected  in  the  range 


69 


70 


FIG  IV-6  RANGE  STATISTICS  FOR  CASE  ONE  OF  EXANPLE  TUO 
SHOUING  THE  ESTINATO'?  WITCHING 


(Q)Hinuizu 


FIG  IU-7  AZinUTH  STATISTICS  FOR  CASE  ONE  OF  EXAMPLE  TUO 


u. 
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■IG  IU-8  ELEUATION  STATISTICS  FOR  CASE  ONE  OF  EXANPLE  TUO 


(QIN0I1UA713 


residuals.  Originally  the  switching  scheme  included  the 
azimuth  and  elevation  residuals  but  the  final  version  uses 
only  the  range  residuals.  For  comparison  purposes,  similar 
plots  for  case  two  and  case  three  are  given  in  Figures  TV-10 
through  IV-15.  The  estimator  switching  algorithm  didn't 
work  as  well  in  either  of  these  cases  as  it  did  in  case  one. 
With  30  points  per  batch,  the  program  switched  from  the  six 
state  estimator  to  the  BRV  estimator  at  185  seconds  then  im¬ 
mediately  switched  to  the  MaRV  estimator.  With  10  points 
per  batch  only  one  switch  is  made  also  this  occurring  at 
192  seconds  when  the  maneuver  is  detected.  These  comparisons 
suggest  the  desirability  of  using  the  larger  batch  sizes. 

Next  a  comparison  of  the  accuracies  of  the  estimated 
parameters  for  the  three  cases  was  made.  This  was  accom¬ 
plished  by  comparing  the  differences  between  the  estimated 
parameters  and  their  true  values.  In  Figures  IV-16  through 
IV-21  the  .’ifference  between  the  estimated  position  and  veloc 
ity  states  and  their  true  values  from  a  single  sample  run 
of  case  one  are  plotted.  In  Figures  IV-22,  IV-23,  arid  IV-24 
the  estimated  values  of  hypersonic  ballistic  coefficient, 
and  are  plotted.  Also  in  each  of  the  Figures  ZV-16 
through  IV-24,  plus  and  minus  3-o  curves  from  the  estimator 
covariances  are  plotted.  In  an  effort  to  nave  space,  when 
the  case  two  and  case  three  comparisons  aro  made  only  the  X 
component  of  position  and  velocity  and  the  aerodynrimic  poro- 
meters  are  included.  Those  are  included  as  riguics  IV-25 
through  IV-34.  The  comparinontt  may  be  nummarlzncl  as  followm 


(a)  Except  for  poaltion,  tha  atata  aatlmataa  ara 
aignlflcantiy  dagrada<i  at  tha  numbar  of  pointa 
par  batch  ia  raduoad  from  50  to  30  to  10. 

(b)  Tha  arrora  in  tha  aatjinataa  of  tha  aarodynamic 
parametara  ara  moat  affactad  by  tha  raduction  in 
tha  numbar  of  pointa  par  batch. 

(c)  In  all  caaaa,  axcapt  juat  baforo  aatimator  awltch- 
ing,  tha  true  error  ia  oontainad  within  the  piua 
and  minua  3a  curvea. 

Notai  tha  apparent  diverganca  of  tha  aatimator  atandard  devia¬ 
tion  on  each  of  tha  Eiguroa  IV-16  through  lV-34  ia  a  raault 
of  the  collapaing  of  tha  batch  aiaa  onca  Uia  laet  batch  of 
data  ia  ancountarad.  For  axampla,  aay  tha  bitch  aiaa  ia 
apacifiad  to  be  50  pointa.  Onoe  tha  laat  50  pointa  of  data 
aro  procaaaed,  than  the  batch  aiaa  la  raduoad  by  tha  numbec 
of  pointa  tha  apooh  ia  advanoad  until  tha  appc)i  raachaa  tha 
final  point.  An  alternate  way  to  maohoniice  thla  would  i>o  to 
juat  taHe  tha  aatlmataa  from  the  laat  full  batch  and  propagate 
them  to  the  final  tim#« 

The  raaulta  of  tha  aingle  aumpla  runa  indioate  that 
tha  atata  paramatera  can  be  aatimated  exceptionally  wall  for 
batch  Nieaa  of  50  pointa.  In  the  region  of  peaK  deceleration! 
where  the  aerodynamic  paranietera  are  estimeted  moat  aocuiitely, 
the  atandard  davl/itlont  aa  determined  from  the  eat  Imator  co- 
variance  matrix  ate  well  within  a  fow  percent  of  the  true 
pacamotor  vrlunM,  ^uoul  2%  for  hdJJlet.lc  cuefflclMiil.  and 
for  Die  turn  and  c/jlmli  parntnot  vi  a .  In  uidni  lo  voiJfy  ihn 
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TO  9-STATE  ESTIHOTOH 


FIC  IV-t9  fMHCE  STftTISTICS  CASE 
S40UINC  THE  ESTI?»(rOll  SUITi 


FIG  lU-ia  ELEUATION  STATISTICS  FOR  CASE  TUO  OF  EXAMPLE  TUO 


riG  IU-13  RANGE  STATISTICS  FOR  CASE  THREE  OF  EXAMPLE  TUO 
SHOUING  THE  ESTIMATOR  SUITCHING 
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FIG  IV-14  AZIMUTH  STATISTICS  FOR  CASE  THREE  OF  EXAMPLE  TUO 


ELEUATION  STATISTICS  FOR  CASE  THREE  OF  EXAMPLE  TUO 


RUE  ERROR  AND  3-SIGNA  ERROR  FOR  A  CASE-ONE 
SINGLE  SANPLE  RUN 


EXARPLE-TUO  SINGLE  SANPLE  RUN 


1 1 


EXAnPLE-TUO  SINGLE  SAMPLE  RUN 


EXAMPLE-TUO  SINGLE  SAHPLE  RUN 


FIG  lU-ae  V  UELOCITV  TRUE  ERROR  AND  S-SIGtlA  ERROR  FOR  A  CASE-ONE 
EXAMPLE-TUO  SINGLE  SAflPLE  RUN 


100.0 
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-22  BALLISTIC  COEFFICIENT  TRUE  ERROR  AND  a-SIGMA  ERROR  FOR  A 

CASE-ONE,  exauple-tuo  single  sample  run 


6.0- 


‘isJI 


PARAMETER  TRUE  ERROR  AND  S-SIGMA  ERROR  FOR 
ONE,  EXAMPLE-TUO  SINGLE  SAMPLE  RUN 


-O’* 


PARAPIETER  TRUE  ERROR  AND  3-SIGnA  ERROR  FOR  A 
ONE,  EXAPIPLE-TUO  SINGLE  SAfIPLE  RUN 


100.0 


EXAMPLE-TUO  SINGLE  SAMPLE  RUN 


100.0 


EXAMPLE-TUO  SINGLE  SAMPLE  RUN 


100.0- 


X  UELOCITV  TRUE  ERROR  AND  3-SIGnA  ERROR  FOR  A  CASE-TUO 
EXAMPLE-TUO  SINGLE  SAMPLE  RUN 


O'OOt 


EXAMPLE-TUO  SINGLE  SANPLE  RUN 


IS.O 


FIG  IU-29  BALLISTIC  COEFFICIENT  TRUE  ERROR  AND  S-SIGNA  ERROR  FOR  A 
CASE-TUO-  EXAPIPLE-TUO  SINGLE  SANPLE  RUN 


>2.0 
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BALLISTIC  COEFFICIENT  TRUE  ERROR  AND  3-SICrJA  ERROR  FOR  A 
CASE-THREE,  EXANPLE-TUO  SINGLE  SANPLE  RUN 
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PARAMETER  TRUE  ERROR  AND  3-SIGnA  ERROR  FOR  A 
TUO,  EXAMPLE-TUO  SINGLE  SAMPLE  RUN 


10.0 


PARAHETER  TRUE  ERROR  AND  S-SIGNA  ERROR  FOR  A 
THREE.  EXAMPLE-TWO  SINGLE  SAMPLE  RUN 


FIG  IU-33  CLIHB  PARAHETER  TRUE  ERROR  AND  S-SICNA  ERROR  FOR  A 
CASE>TUO«  EXARPLE-TUO  SINGLE  SANPLE  RUN 


results  of  the  single  sample  covariance  analysis,  a  ten  run 
Monte  Carlo  simulation  was  performed  for  the  case  one  runs 
(50  point  batch  sizes).  The  results  Of  this  Monte  Carlo 
simulation  are  plotted  in  Figures  IV-35  through  IV-43 .  In 
each  of  these  figures  both  the  standard  deviations  from  the 
Monte  Carlo  simulation  and  from  a  single  sample  run  are 
plotted.  Plus  and  minus  values  are  plotted  so  as  to  locate 
the  zero  line  easily.  The  results  of  this  Monte  Carlo  simula¬ 
tion  show  that,  in  general  the  single  sample  standard  devia¬ 
tions  are  realistic,  albeit  conservative,  estimates  of  the 
errors  in  the  estimated  parameters. 


X  POSITION  COnPARISON  OF  THE  SINGLE  SAM 
DEUIATION  WITH  THE  RMS  ERROR  FROM  A  10-1 
SIMULATION,  EXAMPLE  T(JO,  CASE  ONE  TYPE 
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FIG  IU-36  V  POSITION  CONPARISON  OF  THE  SINGLE  SAMPLE  STANDARD 

DEUIATION  UITH  THE  RMS  ERROR  FROM  A  10>RUN  MONTE  CARLO 
SIMULATION,  EXAMPLE  TUO,  CASE  ONE  TYPE 


2  POSITION  CORPARISON  OF  THE  SINGLE  SANPLE  STANDARD 
DEVIATION  UITH  THE  RNS  ERROR  FRON  A  lO-RUN  NONTE  CARLO 
SIMULATION,  EXAMPLE  TUO,  CASE  ONE  TYPE 


X  UELOCITy  couparison  of  the  single  sanple  standard 
DEVIATION  UITH  THE  RHS  ERROR  FROfI  A  10-RUN  NONTE  CARLO 
SIMULATION,  EXAMPLE  TUO,  CASE  ONE  TYPE 


FIG  IU-39  V  VELOCITY  COnPARISOH  OF  THE  SINGLE  SAMPLE  STANDARD 

DEVIATION  UITH  THE  RI1S  ERROR  FROM  A  IC-RUN  MONTE  CARLO 
SIMULATION,  EXAMPLE  TUO,  CASE  ONE  TYPE 


FIG  IU-40  Z  VELOCITY  COMPARISON  OF  THE  SINGLE  SAMPLE  STANDARD 

DEVIATION  UITH  THE  RMS  ERROR  FROM  A  lO-RUN  MONTE  CARLO 
SIMULATION,  EXAMPLE  TUO,  CASE  ONE  TYPE 


FIG  IU-41  BALLISTIC  COEFFICIENT  CODPARISON  OF  THE  SINGLE  SAMPLE  STANDARD 
DEVIATION  UITH  THE  RMS  ERROR  FROM  A  10-RUN  MONTE  CARLO 
SIMULATION,  EXAMPLE  TUO,  CASE  ONE  TYPE 


4.0 
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FIG  IU-42  TURN  PARAHETER  COMPARISON  OF  THE  SINGLE  SAMPLE  STANDARD 
DEUIATION  UITH  THE  RMS  ERROR  FROM  A  10-RUN  MONTE  CARLO 
SIMULATION,  EXAMPLE  TUO,  CASE  ONE  TYPE 


CLIHB  PARAMETER  COMPARISON  OF  THE  SINGLE  SAMPLE  STANDARD 
DEUIATIOM  WITH  THE  RMS  ERROR  FROM  A  10-RUN  MONTE  CARLO 
SIMULATION.  EXAMPLE  TUO.  CASE  ONE  TYPE 


V  Suininary.  Conclusions  and  Recommendations 


A  very  flexible  algorithm  for  post-flight  analysis  of 
radar. tracking  data  collected  on  a  ballistic  missile  has 
been  presented.  It  was  shown  that  by  specifying  certain 
values  of  input  parameters,  the  algorithm  would  work  in 
either  the  batch  weighted  least  squares  mode  or  as  a  sequen- 
tial  estimator.  Two  numerical  examples  detailing  the  opera¬ 
tion  of  the  estimator  in  both  modes  were  presented.  For  the 
sequential  estimator  mode,  a  method  was  devised  to  allow  the 
program  to  switch  automatically  from  a  free-flight  model  to 
a  BRV  model  and  then  to  a  MaRV  model .  This  means  that  the 
algorithm  can  be  used  for  f roe-flight  or  for  reentry  analyses 
and  the  reentry  may  be  either  BRV  flight  or  MaRV  flight.  A 
particular  state  vector  of  parameters  to  determine  is  selec¬ 
ted  by  initializing  the  state  transition  matrix  to  the  iden¬ 
tity  matrix.  For  each  parameter  of  the  9-element  MaRV  state 
vector  which  is  not  to  be  included  in  the  model,  the  diagonal 
element  of  the  initial  state  transition  matrix  corresponding 
to  that  parameter  is  set  to  zero.  This  allows  much  greater 
flexibility  than  was  actually  demonstrated  in  this  report. 

It  means  that  any  combination  of  the  nine-elements  may  be 
used  as  the  state  vector  of  parameters  tu  be  determined.  As 
an  example,  if  the  position  and  velocity  of  the  RV  were  ac¬ 
curately  known  at  reentry  as  might  be  provided  by  a  long 
free-flight  track,  one  could  conceivably  solve  only  for  the 


aerodynamic  parameters  while  propagating  position  and  velocity 
according  to  the  models  given  in  Chapter  II .  This  could  be 
done  in  either  the  batch  or  sequential' modes .  A  recommenda¬ 
tion  for  future  research  is  to  compare  the  results  of  such  a 
technique  to  those  given  in  this  report. 

The  batch  weighted  least  squares  method  which  was  pre¬ 
sented  in  Example  One  of  Chapter  IV,  offers  such  attractive 
features  when  compared  with  the  sequential  estimator  that 
one  concludes  that  when  the  dynamic  and  measurement  models 
are  accurately  known,  the  batch  method  is  preferred.  The 
attractive  features  mentioned  above  include i 

(a)  running  time  for  the  batch  estimator  is  consider¬ 
ably  less  than  that  required  for  the  sequential 
estimator, 

(b)  the  estimator  uncertainties  are  much  smaller  for 
the  batch  estimator  as  indicated  by  a  comparison  of 
Examples  One  and  Two  of  Chapter  IV,  and 

(c)  programming  of  the  batch  estimator  is  much  simplier 
and  therefore  less  time  consuming  with  few  chances 
of  making  errors. 

These  features  of  the  batch  estimator,  motivate  one  to  look 
^or  ways  of  improving  models,  so  that  the  batch  estimator  is 
applicable.  This  is  probably  not  possible  for  the  MaRV  since 
there  is  no  accurate  general  model  of  maneuvering  flight) 
however,  for  BRV  flight,  it  should  be  possible  to  come  up 
with  some  analytic  model,  perhaps  based  on  wind  tunnel  data, 
to  model  accurately  atmospheric  drag.  One  such  model  was 


suggested  in  Chapter  II  where  was  given  as  a  function  of 
Mach  number 

^  ^  ^  (11-71  ) 

In  the  batch  mode  )<2  and  could  be  included  as  param¬ 

eters  to  be  determined  In  the  estimation  process .  The  drag 
model  given  by  Eg  (11-71)  may  not  be  the  best  model  in  gen¬ 
eral*  A  recommendation  for  further  study  is  to  investigate 
models  like  Eq  (11-71 )  for  the  BRV  case  to  see  under  what 
circumstances  the  batch  estimator  would  give  reasonable  re¬ 
sults.  It  seems  that  this  could  only  be  done  using  real 
data. 

A  final  recommendation  for  future  study  is  in  the  area 
of  sequential  estimation  with  small  batch  sizes.  It  was 
mentioned  earlier  that  for  small  batch  sizes,  a  priori  sta¬ 
tistics  had  to  be  used  to  keep  non-singular  and 

therefore  invertible.  It  would  appear  that  for  very  large 
maneuvers  or  where  and  vary  with  time,  small  batch 
sizes  must  be  considered.  When  they  were  considered  briefly 
in  this  study  it  appeared  that  there  was  an  inverse  relation¬ 
ship  between  the  batch  size  and  a  weighting  factor  applied 
to  the  inverse  of  the  a  priori  covariance,  .  More  re¬ 
search  may  reveal  a  relationship  so  that  the  weighting  can 
be  applied  automatically  by  the  program  from  the  batch  size. 
This  investigation  could  also  look  at  adaptively  optimizing 


and  changing  the  batch  size  as  a  function  of  time.  As  an 
example,  it  may  be  advantageous  to  use  smaller  batch  sizes  in 
regions  of  high  deceleration  and  larger  batch  sizes  in  re¬ 
gions  where  the  deceleration  is  not  as  high. 
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APPENDIX  A  The  Geometrical  Shape  of  the  Earth 

This  appendix  describes  the  model  used  for  the  shape 
of  the  earth  and  developes  the  relationships  n'icessary  in 
seme  of  the  transformations  involved  in  the  trajectory  esti¬ 
mation  problem. 

The  geometrical  shape  of  the  earth  is  modeled  as  an 
ellipsoid  of  revolution  with  the  revolution  being  taken 
about  the  minor  axis  which  coincides  with  the  earth's  axis 
of  rotation.  This  ellipsoid  of  revolution,  which  differs 
only  slightly  from  a  sphere,  is  used  to  locate  positions  on, 
or  above,  the  earth  in  a  geodetic  (or  geographic)  coordinate 
system.  Figure  A-1  is  a  graphic  portrayal  of  the  reference 
ellipsoid  with  the  eccentricity  greatly  exaggerated  showing 
the  parameters  associated  with  the  ellipsoid.  The  values 
given  for  the  semi -major  and  semi -minor  axes  are  from  the 
World  Geodetic  System  (WGS-72)  reference  ellipsoid.  From  the 
geometry  shown  in  figure  A-2,  the  relationship  between  geo¬ 
detic  and  geocentric  latitude  is  determined.  In  the  local 
cartesian  coordinate  system  of  Fig  A-2,  the  equation  of  the 
ellipse  is 


The  geodetic  latitude  line  is  perpendicular  to  the  tangent 
plane  at  the  earth's  surface  as  shown,  so  that 
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a  -  SEMI-MAJOR  AXIS,  6378135(M) 

6  -  SEMI-MINOR  AXIS,  6356750(H) 

A  -  (5E0CENTRIC  ALTITUDE 

C  -  LOCAL  RADIUS  OF  THE  EARTH 

A  -  LONGITUDE 

<t’  -  GEODETIC  LATITUDE 

.<<  -  GEOCENTRIC  LATITUDE 


FIG  A-1  PARAMETERS  ASSOCIATED  UITH  REFERENCE  ELLIPSOID 
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a  -  SEPII-HAJOR  axis,  6378135 (M) 

b-  SEni-niNOR  AXIS,  6356750(n) 

h  -  GEOCENTRIC  ALTITUDE 

«  -  LOCAL  RADIUS  OF  THE  EARTH 

0  -  LONGITUDE 

f’-  GEODETIC  LATITUDE 

4’-  GEOCENTRIC  LATITUDE 


FIG  A-a  CROSS  SECTION  OF  THE  REFERENCE  ELLIPSOID 
USED  FOR  DETERMINING  THE  RELATIONSHIP  BETUEEN 
GEODETIC  AND  GEOCENTRIC  LATITUDE 
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(A-2) 


Differentiating  equation  (A-1)  we  get 


(A-3) 


Substituting  (A-2)  in  (A-3),  ve  obtain 


(A-4) 


From  the  geometry  of  figure  (A-2),  we  have 


TAf^  <9 


(A-5) 


Putting  (A-5)  into  (A-4)  we  obtain 


(A-6) 


which  is  the  relationship  between  geodetic  and  geocentric 
latitude  which  was  sought. 

The  radius  of  the  earth  at  CJ  is  also  determined  from 
the  geometry  of  Fig  A-2. 


4  X 

^  =  X 


(A-7) 
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Solve  (A-1)  for  y  and  substitute  in  (A-7 )  to  obtain 


Divide 


through  by 


2 


to  obtain 


Rearranging  (A-9)  and  recognizing 


CQS^(^ 


we  obtain 


(A-8) 


(A-9) 


(A-10) 


which  gives  the  local  radius  of  the  earth  for  any  geocentric 
latitude  where  e  is  the  eccentricity  given  by 

(A-12) 
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APPENDIX  B  Partials  Required  in  the  Differential 

Correction  Process 

Central  to  the  working  of  the  differential  correction 
process  are  the  partial  derivatives  of  the  observations  -  in 
this  case  range  (R),  azimuth  (A),  and  elevation  (E)  -  with 
respect  to  the  parameters  to  be  estimated,  in  this  case  the 
elements  of  the  state  vector  at  some  epoch  time,  t^.  The 
most  general  state  vector  for  a  MaRV  is  nine  dimensional. 

Let  the  state  vector  at  the  epoch  time  be  defined  as 

Cr,Cc)^  (B-1) 

where  the  nine  elements  are  as  previously  defined.  No  o  sub¬ 
script  is  required  for  and  C^  since  they  are  modeled  as 
constants.  The  partials  required  for  the  differential  cor¬ 
rection  process  then  may  be  symbolized  as 


dR. 
_ 1* 

tfi. 

3r. 

_ If 

aR. 

_ If 

dR. 

_ If 

dR. 
_ 1 

a  c 

^T 

^Cc 

dA. 
_ 1* 

aA. 

_ If 

aA. 

_ If 

aA. 

_ If 

aA . 

_ If 

aA. 

_ If 

aA. 

_ If 

dA. 
_ 1 

'^o 

3Z 

o 

o 

a  c 

T 

ifi. 

aE. 

_ If 

as, 

_ If 

aE. 

_ If 

as. 

_ If 

aE. 

_ If 

aE. 

_ If 

dE. 
_ 1 

ax 

o 

ay 

o 

^'=T 
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where  the  subscript  i  is  used  to  denote  the  time,  t^  of  the 
i^^  measurement  i.e.,  =R(t^).  If  there  are  n  measure¬ 

ments,  then  i  =  1,  2,  3,  . . .  n.  tj^  may  or  may  not  equal  t^. 
In  this  developme  t  t^^  is  always  greater  than  or  equal  t^. 

The  radar  range,  azimuth,  and  elevation  at  t^  are  all 
functions  of  the  RV’s  position  at  t^^. 


(B-2) 

4  "  V/,  Zi ) 

(B-3) 

fiC A,  Vi,!,) 

(B-4) 

where  X^,  Y^,  and  are  elements  of  the  state  a  t^^  in  the 
ECI  coordinate  system, 

(B-2)  is  determined  from  (B-3)  by  the  chain  rule  for  partial 
derivatives  given  by  the  following  matrix  equation 


'hi’hu 

Wv/,  u/i/* 

—  — 

'ifi/iZ, 

ziihi* 

Wik,  ^Z/M 

'bii/iL 

'iKi/Xf 

'hii/lt  'his /lit* 

b/,-A/4  az,'/V. 

iXi/ACc  WvACc  aZj/aCc 

where  the  matrix  is  a  partitian  of  the  state  transition 
matrix  )  to  be  developed  on  the  next  few  pages. 

Azimuth  and  elevation  partials  are  calculated  similar 
to  (B-4)  using  the  same  partition  of 
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The  State  Transition  Matrix 

The  state  transition  matrix  is  defined  as 


^ 

VC,  JX,  9X,  J/, 

^  ij  Sli  ^  M 

n  >4  H  Tx, 

Ui  i/i  Si  M  jii  :& 

ih,  32.  32.  32,  Vi  32, 

wc  >>;•  «i  ii  jy,  ii- 

3X,  3X,  3X.  >X, 

„  ,  M  &  &  J  a 

=  if.  n.  )K  w.  K  », 

a  wi-  S-  %  4  & 

JZo  3Z*  ^  32,  Tz.  32, 

iJf.’  iX-  ii-  jX;  3l 

}A  3/« 

Ml  M  IZi  >/  iJ-  * 

3^  3^/"  3^  2Cr  3^r 

JDi  32/  3^  3^- 

3Ce  3Cf  3<t  3^ 


^  ^  3k 
>/,  “BX,  3X. 

^  jr, 

tx  3/,  3/. 

3Z.  32,  "32, 

4^.- 

T3r.  3x.  3Xe 
Ty.  3X 

3i.  32.  32, 

3^  ^  3A 

3<rr  3t^  Kr 
^4  s<:c 


The  state  transition  is  propagated  to  the  time  of  each  ob¬ 
servation  by  the  matrix  chain  rule  equation 


A.)  '^0/,  i.  ) 


(B-7) 


^^4  -T 


(B-8) 


where  ^  is  the  identity  matrix  and  the  elements  of 

* 

are  determined  from  the  truncated  Taylor's  series 


Z 


(B-9) 


These  simplifying  approximations  are  used  only  for 
partial  calculations  and  not  for  propagating  the  state. 
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i- 


P- 

x,v, 

(B-10) 

1 

f-- 

i 

r 

i 

(B-11) 
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t 

•  •  ••  ■ 

Xv, 

(B-12) 

ir. 

M- 

f 

Zi„-  i-i'ziii' 

(B-13) 

t 

(B-14) 

P 

# 

f; 

t 

(B-15) 

II 

(B-16) 

V  vl- 

r 

& 

(B-17) 

w. 

•  «  •  «  •  • 

f 

where  X,  Y,  and  Z  are  the  components  of  the  total 

accelera- 

tion  due  to  gravity,  drag,  turn,  and  climb  given  in  Chc..-)ter 


Partials  of  Gravity  Acceleration 

The  earth's  gravity  model  which  is  used  in  propagahing 
the  state  vector  is  given  in  Chapter  II }  however,  a  simpler 
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approximation  is  sufficiently  accuratfe  for  calculating  the 
partials 


(B-1 

(B-2 

>  r 

(B-2 

where 

/S' 

(B-2 

Then 

>V«- 

(B-2 

^ 

(B-2 

^^/dZi  * 

(B-2 

4/j)0  -  ’Vs/ 

( B— 2 

(B-29) 


=  (B-30) 

All  other  partials  of  the  gravity  acceleration  with  respect 
to  the  state  are  2ero. 

Partials  of  Drag  Acceleration 

The  partials  of  acceleration  due  to  drag  are  determined 
from  the  following  equations 


(B-32) 

•*  **  V  ' 

(B-33) 

(B-34) 

(B-35) 

y-i= 

(B-36) 

(B-3  7) 
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(B-38) 


The  twenty-seven  partials  of  drag  acceleration  with  respect 
to  the  state  elements  are  then 


(B-39) 

(B-40) 

3X/  ax/  )u,{ 

(B-41) 

a/t  V-  ^  ^Wo\r^ 

(B-42) 

(B-43) 

av,  av,.  y,v 

(B-44) 

^  4^-  i^' 

aZ/  ^  1/5,  7000/; 

(B-45) 

?Zt  32.|  Xl*^' 

(B-46) 

§1-/, 

3Z*  4^' 

(B-47) 

lli  ait  t,-  '^'1^ 

(B-48) 
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(B-49) 

(B-50) 

ay,  u.' 

(B-51) 

(B-52) 

3^/.  2LL 

(B-53) 

^4^/ 

92/  2X/ 

(B-54) 

II 

(B-55) 

2^/— 

3Zi'  3^1' 

(B-56) 

3>k  r  - 
?A  A 

{B-57) 

3^-*  - 
3/3L»  /»* 

(B-58) 

?A  A 

(B-59) 

Partials  of  Turn  Acceleration 

The  components  of  the  acceleration  due  to  turn  were 
developed  in  Chapter  II  as 


Tt 

(B-60) 

%  ’CrA  VilMk  +-^X)  -Xitii 

1 

(B-61) 

Zr.’-CrA  -%!X  f/lX)] 

7 

(B-62) 

(B-63) 

In  an  effort  to  simplify  the  notation  in  the  rather  complicated 
turn  partials,  the  i  subscript  vhich  denotes  those  parameters 
which  vary  with  time,  will  be  dropped.  Then  the  partials  of 
turn  acceleration  with  respect  to  the  nine  state  elements 
are 


{B-64) 

(B-65) 

(B-66) 

(B-67) 

(B-68) 
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(B-69) 


Hi  T  r‘  J  T  31 


(B-72) 

(B-73) 

(B-74) 

(B-75) 

(B-76) 

(B-77) 

(B-78) 

(B-79) 

(B-80) 


(B-81) 

Where 

(yi-zVza  -(z^-)(z)z 

2x  r 

(B-82) 

ir»  (yz-zi)z  iCzjcrxzJzn  -r/i-yx)/'>L^-n-y) 

(B-83) 

(B-84) 

jy  r 

ar  «  -(yi-2y)t  ^  fzirxz)x^ 

T 

2z:  =  (z^-A'z)z -Cx)i.-vo<0/ 

r- 

(B-85) 

iT -  “(yz-i^Z’^ LtXry^x 

jy  7- 

(B-86) 

21  _  ^v2.-7^)y-(z;t-xz)x 

3Z.  7-  ■  ■ 

(B-87) 

and  where 

^  ^  ) 

(B-88) 

(B-89) 

.  -/D^vC 

^  2,  ^000  r 

(B-90) 
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(B-91) 


9 


s  ;8^Z 


(B-92) 


(B-93) 


Additionally 

3Xr  - 

TtCj  '  r 

(B-94) 

3Yr  _  ^K’(zat.-xz) 

acr  ^ 

{B-95) 

iiz 

liCf  -r 

(B-96) 

•  •  « •  •• 

Partials  of  X^,  and  wrt  8^  and  are  all 

zero. 

Partials  of  Climb  Acceleration 

The  components  of  the  acceleration  due  to  climb  from 
Chapter  II  are 


(B-97) 


(B-98) 
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^)'Ya,CY2'ZyJ 


from  which  the  partials  of  the  climb  acceleration  with  respect 
to  the  state  elements  are 

.  %. - Y>^)  ~ Z C^X*r yz )  (B-lOO) 

T 


101) 


3i  -^ii(zi*xk.)]'[zCyz-zt)  -)^CxYt-Y;(.)]|fj 

^  iC/i  I  (B-ioi) 

,  X*Cz^'’XZ)'’^^y2 -zi.)  1  (b-102) 

T  -^x  j 

^  -c  ’Z(z^-<xz3^J 

f  -yx>)  -zCzx^~xz) 


(B-103) 


-z(z^'X 


MS) 
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2( yz -/iu.)  - -  yx«)  3(g(4) 
T 


Z’*Sll22X^-Xi)j-\jU.CzX^-XZ) 


LOr)(^-xz) -tCvz-zt')  d(M 
T  "ay 


^  -Li(xi-Y^)  -z(zx^-xz)]^j 

^  t(xy4.-y;l)-iCzywz)  3(^|^)| 
-riiiZ  -  [z  (y'z  -Siz)  -  v^)]  IF 


3Z 
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2X1 

21 


+  z( yi  -  zV)  (xj  -  yp  9  (^1^  \ 

r  '^^  j 


^  ^  L  \r(xl^ </:)  -[Mzx.- xz)  -  y^vz  'ZY.^Ifj 


^  >u(zi.-xz)-t(y^-zt)  dC>ti) 
- r - 


ii 

3;i: 


I .  r.  U 1^  r^yXu-  Zi)  -[>:  ft  V  v^)  -  i  |f  j 


r 


(B-104) 


(B-105) 


(B-106) 


(B-107) 


(B-108) 


'ii(xt'yk^)-zCzt-XZ)  d(^l^]\  (b_109) 


I  S 

3X~  j 
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^.z(yz-z9.}->:,(xt-yL)  a^Jl  (b-uo) 

r  3>!  J 

j^7-(2Z)^-xz)  -  L^(zx^-/z)  -?^(yi-zyo]^J 

^  ^(zL'Xz)-X(yi-z9^  2^l^)|  (b-hd 
;  -  y^)  -  K  (xl-  y^)  -  i  -X  z  j 


(B-112) 


^  -[2cyz  -zi)  -;Lcy)l-yA.)]|fj 


,  z(yi-z9)-L(.x9^-y)0  (B-113) 


zzy:)-[_^iz'>^-xt)-t 


Uvz  'ZtiU' 


-!i^(yz~z^)  (b-114) 


c.Up:lk 


-^xz)-r'2(xl-yx»)- 
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X>)-Z(ZVX1^' 


)i.(xt-yL')''tCz:u-xl)  ^ 


(B-115) 


^  t  •■ - :  tr^  «-  —1^%, 


^ . Q U  'ruyz-zi.)  -[i 


-lz(vz-zt)  '^(xi-yL)] 


T  ^ 


(B-116) 


^ -(^^i/r''(-xx:-Y»)-  r  Uz^-xz)  (yz-z)i}|?  j 


jL(z^-xz)-y^(y2-2i:) 


(B-117) 


where  the  partials  of  T  are  as  given  in  the  previous  section 


•ax  nooor 

'w  'ir  •jooof 

37  '  7<S>0<7/^ 


'h(fp  .  p  X. 


(B-118) 


(B-119) 


(B-120) 


(B-121) 


(B-122) 


{L-123) 


The  partials  of  ,  and  wrt  8^  and  are  all  zero 
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d ..  aii:  FMiiMjiiifMi}' 
Jt  L  r 


(B-124) 


?^-^t:rz(Yi-7i0-^(x-fc-v>O 
ja  L  T 


(B-125) 


(B-126) 


With  the  partials  of  acceleration  with  respect  to  state  ele¬ 
ments  given  in  this  and  previous  sections  and  equations  (B-9) 
through  (B-18)  all  of  the  elements  of  hand. 
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neuvering  reentry  vehicle  (MaRV)  from  post-flight  analysis  of  rad 
tracking  data  was  developed.  The  technique  uses  a  sequential 
weiglited  least  squares  estimator  which  steps  through  the  data 
processing  batches  of  data  sequentially.  The  observations 
consisted  of  angle  and  range  measurements  from  a  precision 
tracking  radar  located  in  the  vicinity  of  the  impact  area.  -  ^ 
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A  SIX  dimension  estimator  consisting  of  three  position  and 

three  velocity  components  was  designed  for  estimating  the 
state  of  the  RV  during  free  flight.  Two  reentry  estimators 
were  developed.  One,  a  seven  dimension  estimator  consisting 
of  the  six  elements  of  the  free-f light  estimator  plus  the 
ballistic  coefficient,  was  designed  to  estimate  the  state 
during  ballistic  reentry.  The  second  is  the  nine  element 
MaRV  estimator  used  once  a  maneuver  is  begun.  An  algorithm 
based  on  measurement  residual  monitoring  was  developed  to 
switch  adaptively  from  the  six-state  estimator  to  the  seven- 
state  estimator  and  then  to  the  nine-state  estimator.  A 
series  of  numerical  simulations  was  performed  to  validate 
the  technique  and  its  programming.  Monte  Carlo  simulations 
were  used  to  verify  the  accuracy  of  the  estimator  covariance 
matrix . 
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